tyleransom/autochoice.csv

## autochoice.csv

          
            id
            male
            income
            size
            choicecat
            dealer1
            dealer2
            dealer3

            
              1
              1
              46.7
              2
              3
              18
              8
              5

            
              2
              1
              26.1
              3
              1
              17
              6
              2

            
              3
              1
              32.7
              4
              1
              12
              6
              2

            
              4
              0
              49.2
              2
              2
              18
              7
              4

            
              5
              1
              24.3
              4
              1
              10
              7
              3

            
              6
              0
              39
              4
              1
              13
              4
              1

            
              7
              1
              33
              2
              1
              24
              7
              4

            
              8
              1
              20.3
              4
              1
              9
              3
              3

            
              9
              1
              38
              3
              2
              14
              10
              2

            
              10
              0
              60.4
              3
              1
              18
              5
              4

            
              11
              1
              69
              3
              1
              23
              8
              2

            
              12
              1
              27.7
              1
              2
              21
              9
              5

            
              13
              1
              41
              2
              2
              20
              10
              5

            
              14
              1
              65.6
              4
              2
              13
              4
              5

            
              15
              1
              24.8
              3
              1
              21
              5
              3

            
              16
              1
              21.6
              2
              1
              19
              8
              3

            
              17
              0
              43.3
              1
              1
              23
              9
              5

            
              18
              0
              44.4
              1
              2
              22
              10
              4

            
              19
              1
              45.3
              4
              1
              17
              8
              1

            
              20
              1
              48.8
              3
              1
              20
              10
              5

            
              21
              1
              44.7
              4
              3
              23
              7
              2

            
              22
              1
              57.5
              1
              1
              24
              9
              4

            
              23
              1
              42
              2
              1
              20
              8
              3

            
              24
              1
              40.5
              1
              1
              23
              10
              4

            
              25
              1
              44.4
              1
              1
              24
              10
              5

            
              26
              0
              50.3
              2
              1
              21
              7
              4

            
              27
              0
              41.1
              4
              1
              13
              3
              2

            
              28
              1
              49
              3
              2
              22
              10
              3

            
              29
              0
              43
              3
              1
              21
              7
              2

            
              30
              1
              43.8
              4
              3
              16
              8
              3

            
              31
              1
              46.6
              2
              3
              24
              8
              5

            
              32
              1
              45.7
              3
              1
              14
              10
              2

            
              33
              1
              69.8
              3
              1
              23
              10
              5

            
              34
              0
              45.6
              2
              2
              21
              7
              5

            
              35
              1
              20.9
              1
              1
              23
              10
              4

            
              36
              1
              30.6
              4
              1
              20
              4
              4

            
              37
              0
              40.7
              2
              3
              18
              7
              3

            
              38
              1
              46.7
              2
              3
              18
              8
              5

            
              39
              1
              26.1
              3
              1
              17
              6
              2

            
              40
              1
              32.7
              4
              1
              12
              6
              2

            
              41
              0
              49.2
              2
              2
              18
              7
              4

            
              42
              1
              24.3
              4
              1
              10
              7
              3

            
              43
              0
              39
              4
              1
              13
              4
              1

            
              44
              1
              33
              2
              1
              24
              7
              4

            
              45
              1
              20.3
              4
              1
              9
              3
              3

            
              46
              1
              38
              3
              2
              14
              10
              2

            
              47
              0
              60.4
              3
              1
              18
              5
              4

            
              48
              1
              69
              3
              1
              23
              8
              2

            
              49
              1
              27.7
              1
              2
              21
              9
              5

            
              50
              1
              41
              2
              2
              20
              10
              5

            
              51
              1
              65.6
              4
              2
              13
              4
              5

            
              52
              1
              24.8
              3
              1
              21
              5
              3

            
              53
              1
              21.6
              2
              1
              19
              8
              3

            
              54
              0
              43.3
              1
              1
              23
              9
              5

            
              55
              0
              44.4
              1
              2
              22
              10
              4

            
              56
              1
              45.3
              4
              1
              17
              8
              1

            
              57
              1
              48.8
              3
              1
              20
              10
              5

            
              58
              1
              44.7
              4
              3
              23
              7
              2

            
              59
              1
              57.5
              1
              1
              24
              9
              4

            
              60
              1
              42
              2
              1
              20
              8
              3

            
              61
              1
              40.5
              1
              1
              23
              10
              4

            
              62
              1
              44.4
              1
              1
              24
              10
              5

            
              63
              0
              50.3
              2
              1
              21
              7
              4

            
              64
              0
              41.1
              4
              1
              13
              3
              2

            
              65
              1
              49
              3
              2
              22
              10
              3

            
              66
              0
              43
              3
              1
              21
              7
              2

            
              67
              1
              43.8
              4
              3
              16
              8
              3

            
              68
              1
              46.6
              2
              3
              24
              8
              5

            
              69
              1
              45.7
              3
              1
              14
              10
              2

            
              70
              1
              69.8
              3
              1
              23
              10
              5

            
              71
              0
              45.6
              2
              2
              21
              7
              5

            
              72
              1
              20.9
              1
              1
              23
              10
              4

            
              73
              1
              30.6
              4
              1
              20
              4
              4

            
              74
              0
              40.7
              2
              3
              18
              7
              3

            
              75
              1
              46.7
              2
              3
              18
              8
              5

            
              76
              1
              26.1
              3
              1
              17
              6
              2

            
              77
              1
              32.7
              4
              1
              12
              6
              2

            
              78
              0
              49.2
              2
              2
              18
              7
              4

            
              79
              1
              24.3
              4
              1
              10
              7
              3

            
              80
              0
              39
              4
              1
              13
              4
              1

            
              81
              1
              33
              2
              1
              24
              7
              4

            
              82
              1
              20.3
              4
              1
              9
              3
              3

            
              83
              1
              38
              3
              2
              14
              10
              2

            
              84
              0
              60.4
              3
              1
              18
              5
              4

            
              85
              1
              69
              3
              1
              23
              8
              2

            
              86
              1
              27.7
              1
              2
              21
              9
              5

            
              87
              1
              41
              2
              2
              20
              10
              5

            
              88
              1
              65.6
              4
              2
              13
              4
              5

            
              89
              1
              24.8
              3
              1
              21
              5
              3

            
              90
              1
              21.6
              2
              1
              19
              8
              3

            
              91
              0
              43.3
              1
              1
              23
              9
              5

            
              92
              0
              44.4
              1
              2
              22
              10
              4

            
              93
              1
              45.3
              4
              1
              17
              8
              1

            
              94
              1
              48.8
              3
              1
              20
              10
              5

            
              95
              1
              44.7
              4
              3
              23
              7
              2

            
              96
              1
              57.5
              1
              1
              24
              9
              4

            
              97
              1
              42
              2
              1
              20
              8
              3

            
              98
              1
              40.5
              1
              1
              23
              10
              4

            
              99
              1
              44.4
              1
              1
              24
              10
              5

            
              100
              0
              50.3
              2
              1
              21
              7
              4

            
              101
              0
              41.1
              4
              1
              13
              3
              2

            
              102
              1
              49
              3
              2
              22
              10
              3

            
              103
              0
              43
              3
              1
              21
              7
              2

            
              104
              1
              43.8
              4
              3
              16
              8
              3

            
              105
              1
              46.6
              2
              3
              24
              8
              5

            
              106
              1
              45.7
              3
              1
              14
              10
              2

            
              107
              1
              69.8
              3
              1
              23
              10
              5

            
              108
              0
              45.6
              2
              2
              21
              7
              5

            
              109
              1
              20.9
              1
              1
              23
              10
              4

            
              110
              1
              30.6
              4
              1
              20
              4
              4

            
              111
              0
              40.7
              2
              3
              18
              7
              3

            
              112
              1
              46.7
              2
              3
              18
              8
              5

            
              113
              1
              26.1
              3
              1
              17
              6
              2

            
              114
              1
              32.7
              4
              1
              12
              6
              2

            
              115
              0
              49.2
              2
              2
              18
              7
              4

            
              116
              1
              24.3
              4
              1
              10
              7
              3

            
              117
              0
              39
              4
              1
              13
              4
              1

            
              118
              1
              33
              2
              1
              24
              7
              4

            
              119
              1
              20.3
              4
              1
              9
              3
              3

            
              120
              1
              38
              3
              2
              14
              10
              2

            
              121
              0
              60.4
              3
              1
              18
              5
              4

            
              122
              1
              69
              3
              1
              23
              8
              2

            
              123
              1
              27.7
              1
              2
              21
              9
              5

            
              124
              1
              41
              2
              2
              20
              10
              5

            
              125
              1
              65.6
              4
              2
              13
              4
              5

            
              126
              1
              24.8
              3
              1
              21
              5
              3

            
              127
              1
              21.6
              2
              1
              19
              8
              3

            
              128
              0
              43.3
              1
              1
              23
              9
              5

            
              129
              0
              44.4
              1
              2
              22
              10
              4

            
              130
              1
              45.3
              4
              1
              17
              8
              1

            
              131
              1
              48.8
              3
              1
              20
              10
              5

            
              132
              1
              44.7
              4
              3
              23
              7
              2

            
              133
              1
              57.5
              1
              1
              24
              9
              4

            
              134
              1
              42
              2
              1
              20
              8
              3

            
              135
              1
              40.5
              1
              1
              23
              10
              4

            
              136
              1
              44.4
              1
              1
              24
              10
              5

            
              137
              0
              50.3
              2
              1
              21
              7
              4

            
              138
              0
              41.1
              4
              1
              13
              3
              2

            
              139
              1
              49
              3
              2
              22
              10
              3

            
              140
              0
              43
              3
              1
              21
              7
              2

            
              141
              1
              43.8
              4
              3
              16
              8
              3

            
              142
              1
              46.6
              2
              3
              24
              8
              5

            
              143
              1
              45.7
              3
              1
              14
              10
              2

            
              144
              1
              69.8
              3
              1
              23
              10
              5

            
              145
              0
              45.6
              2
              2
              21
              7
              5

            
              146
              1
              20.9
              1
              1
              23
              10
              4

            
              147
              1
              30.6
              4
              1
              20
              4
              4

            
              148
              0
              40.7
              2
              3
              18
              7
              3

            
              149
              1
              46.7
              2
              3
              18
              8
              5

            
              150
              1
              26.1
              3
              1
              17
              6
              2

            
              151
              1
              32.7
              4
              1
              12
              6
              2

            
              152
              0
              49.2
              2
              2
              18
              7
              4

            
              153
              1
              24.3
              4
              1
              10
              7
              3

            
              154
              0
              39
              4
              1
              13
              4
              1

            
              155
              1
              33
              2
              1
              24
              7
              4

            
              156
              1
              20.3
              4
              1
              9
              3
              3

            
              157
              1
              38
              3
              2
              14
              10
              2

            
              158
              0
              60.4
              3
              1
              18
              5
              4

            
              159
              1
              69
              3
              1
              23
              8
              2

            
              160
              1
              27.7
              1
              2
              21
              9
              5

            
              161
              1
              41
              2
              2
              20
              10
              5

            
              162
              1
              65.6
              4
              2
              13
              4
              5

            
              163
              1
              24.8
              3
              1
              21
              5
              3

            
              164
              1
              21.6
              2
              1
              19
              8
              3

            
              165
              0
              43.3
              1
              1
              23
              9
              5

            
              166
              0
              44.4
              1
              2
              22
              10
              4

            
              167
              1
              45.3
              4
              1
              17
              8
              1

            
              168
              1
              48.8
              3
              1
              20
              10
              5

            
              169
              1
              44.7
              4
              3
              23
              7
              2

            
              170
              1
              57.5
              1
              1
              24
              9
              4

            
              171
              1
              42
              2
              1
              20
              8
              3

            
              172
              1
              40.5
              1
              1
              23
              10
              4

            
              173
              1
              44.4
              1
              1
              24
              10
              5

            
              174
              0
              50.3
              2
              1
              21
              7
              4

            
              175
              0
              41.1
              4
              1
              13
              3
              2

            
              176
              1
              49
              3
              2
              22
              10
              3

            
              177
              0
              43
              3
              1
              21
              7
              2

            
              178
              1
              43.8
              4
              3
              16
              8
              3

            
              179
              1
              46.6
              2
              3
              24
              8
              5

            
              180
              1
              45.7
              3
              1
              14
              10
              2

            
              181
              1
              69.8
              3
              1
              23
              10
              5

            
              182
              0
              45.6
              2
              2
              21
              7
              5

            
              183
              1
              20.9
              1
              1
              23
              10
              4

            
              184
              1
              30.6
              4
              1
              20
              4
              4

            
              185
              0
              40.7
              2
              3
              18
              7
              3

            
              186
              1
              46.7
              2
              3
              18
              8
              5

            
              187
              1
              26.1
              3
              1
              17
              6
              2

            
              188
              1
              32.7
              4
              1
              12
              6
              2

            
              189
              0
              49.2
              2
              2
              18
              7
              4

            
              190
              1
              24.3
              4
              1
              10
              7
              3

            
              191
              0
              39
              4
              1
              13
              4
              1

            
              192
              1
              33
              2
              1
              24
              7
              4

            
              193
              1
              20.3
              4
              1
              9
              3
              3

            
              194
              1
              38
              3
              2
              14
              10
              2

            
              195
              0
              60.4
              3
              1
              18
              5
              4

            
              196
              1
              69
              3
              1
              23
              8
              2

            
              197
              1
              27.7
              1
              2
              21
              9
              5

            
              198
              1
              41
              2
              2
              20
              10
              5

            
              199
              1
              65.6
              4
              2
              13
              4
              5

            
              200
              1
              24.8
              3
              1
              21
              5
              3

            
              201
              1
              21.6
              2
              1
              19
              8
              3

            
              202
              0
              43.3
              1
              1
              23
              9
              5

            
              203
              0
              44.4
              1
              2
              22
              10
              4

            
              204
              1
              45.3
              4
              1
              17
              8
              1

            
              205
              1
              48.8
              3
              1
              20
              10
              5

            
              206
              1
              44.7
              4
              3
              23
              7
              2

            
              207
              1
              57.5
              1
              1
              24
              9
              4

            
              208
              1
              42
              2
              1
              20
              8
              3

            
              209
              1
              40.5
              1
              1
              23
              10
              4

            
              210
              1
              44.4
              1
              1
              24
              10
              5

            
              211
              0
              50.3
              2
              1
              21
              7
              4

            
              212
              0
              41.1
              4
              1
              13
              3
              2

            
              213
              1
              49
              3
              2
              22
              10
              3

            
              214
              0
              43
              3
              1
              21
              7
              2

            
              215
              1
              43.8
              4
              3
              16
              8
              3

            
              216
              1
              46.6
              2
              3
              24
              8
              5

            
              217
              1
              45.7
              3
              1
              14
              10
              2

            
              218
              1
              69.8
              3
              1
              23
              10
              5

            
              219
              0
              45.6
              2
              2
              21
              7
              5

            
              220
              1
              20.9
              1
              1
              23
              10
              4

            
              221
              1
              30.6
              4
              1
              20
              4
              4

            
              222
              0
              40.7
              2
              3
              18
              7
              3

            
              223
              1
              46.7
              2
              3
              18
              8
              5

            
              224
              1
              26.1
              3
              1
              17
              6
              2

            
              225
              1
              32.7
              4
              1
              12
              6
              2

            
              226
              0
              49.2
              2
              2
              18
              7
              4

            
              227
              1
              24.3
              4
              1
              10
              7
              3

            
              228
              0
              39
              4
              1
              13
              4
              1

            
              229
              1
              33
              2
              1
              24
              7
              4

            
              230
              1
              20.3
              4
              1
              9
              3
              3

            
              231
              1
              38
              3
              2
              14
              10
              2

            
              232
              0
              60.4
              3
              1
              18
              5
              4

            
              233
              1
              69
              3
              1
              23
              8
              2

            
              234
              1
              27.7
              1
              2
              21
              9
              5

            
              235
              1
              41
              2
              2
              20
              10
              5

            
              236
              1
              65.6
              4
              2
              13
              4
              5

            
              237
              1
              24.8
              3
              1
              21
              5
              3

            
              238
              1
              21.6
              2
              1
              19
              8
              3

            
              239
              0
              43.3
              1
              1
              23
              9
              5

            
              240
              0
              44.4
              1
              2
              22
              10
              4

            
              241
              1
              45.3
              4
              1
              17
              8
              1

            
              242
              1
              48.8
              3
              1
              20
              10
              5

            
              243
              1
              44.7
              4
              3
              23
              7
              2

            
              244
              1
              57.5
              1
              1
              24
              9
              4

            
              245
              1
              42
              2
              1
              20
              8
              3

            
              246
              1
              40.5
              1
              1
              23
              10
              4

            
              247
              1
              44.4
              1
              1
              24
              10
              5

            
              248
              0
              50.3
              2
              1
              21
              7
              4

            
              249
              0
              41.1
              4
              1
              13
              3
              2

            
              250
              1
              49
              3
              2
              22
              10
              3

            
              251
              0
              43
              3
              1
              21
              7
              2

            
              252
              1
              43.8
              4
              3
              16
              8
              3

            
              253
              1
              46.6
              2
              3
              24
              8
              5

            
              254
              1
              45.7
              3
              1
              14
              10
              2

            
              255
              1
              69.8
              3
              1
              23
              10
              5

            
              256
              0
              45.6
              2
              2
              21
              7
              5

            
              257
              1
              20.9
              1
              1
              23
              10
              4

            
              258
              1
              30.6
              4
              1
              20
              4
              4

            
              259
              0
              40.7
              2
              3
              18
              7
              3

            
              260
              1
              46.7
              2
              3
              18
              8
              5

            
              261
              1
              26.1
              3
              1
              17
              6
              2

            
              262
              1
              32.7
              4
              1
              12
              6
              2

            
              263
              0
              49.2
              2
              2
              18
              7
              4

            
              264
              1
              24.3
              4
              1
              10
              7
              3

            
              265
              0
              39
              4
              1
              13
              4
              1

            
              266
              1
              33
              2
              1
              24
              7
              4

            
              267
              1
              20.3
              4
              1
              9
              3
              3

            
              268
              1
              38
              3
              2
              14
              10
              2

            
              269
              0
              60.4
              3
              1
              18
              5
              4

            
              270
              1
              69
              3
              1
              23
              8
              2

            
              271
              1
              27.7
              1
              2
              21
              9
              5

            
              272
              1
              41
              2
              2
              20
              10
              5

            
              273
              1
              65.6
              4
              2
              13
              4
              5

            
              274
              1
              24.8
              3
              1
              21
              5
              3

            
              275
              1
              21.6
              2
              1
              19
              8
              3

            
              276
              0
              43.3
              1
              1
              23
              9
              5

            
              277
              0
              44.4
              1
              2
              22
              10
              4

            
              278
              1
              45.3
              4
              1
              17
              8
              1

            
              279
              1
              48.8
              3
              1
              20
              10
              5

            
              280
              1
              44.7
              4
              3
              23
              7
              2

            
              281
              1
              57.5
              1
              1
              24
              9
              4

            
              282
              1
              42
              2
              1
              20
              8
              3

            
              283
              1
              40.5
              1
              1
              23
              10
              4

            
              284
              1
              44.4
              1
              1
              24
              10
              5

            
              285
              0
              50.3
              2
              1
              21
              7
              4

            
              286
              0
              41.1
              4
              1
              13
              3
              2

            
              287
              1
              49
              3
              2
              22
              10
              3

            
              288
              0
              43
              3
              1
              21
              7
              2

            
              289
              1
              43.8
              4
              3
              16
              8
              3

            
              290
              1
              46.6
              2
              3
              24
              8
              5

            
              291
              1
              45.7
              3
              1
              14
              10
              2

            
              292
              1
              69.8
              3
              1
              23
              10
              5

            
              293
              0
              45.6
              2
              2
              21
              7
              5

            
              294
              1
              20.9
              1
              1
              23
              10
              4

            
              295
              1
              30.6
              4
              1
              20
              4
              4

## test_autochoice_mew1.jl
using DelimitedFiles
using Base.Iterators: flatten
using LinearAlgebra: diag
using Optim
using DataFrames
using HTTP
using uCSV
using LineSearches

# Read in auto data from Stata example
fname = "https://gist.githubusercontent.com/tyleransom/696f3ef806fa4c84e1902f34b8955ca4/raw/49d3f51abe9f126e48907522a69e14b2fd7b5498/autochoice.csv"
data = DataFrame(uCSV.read(IOBuffer(HTTP.get(fname).body), quotes='"', header=1))
data = Matrix(data)
id = data[:,1]
male = data[:,2]
income = data[:,3]
hhsize = data[:,4]
X = hcat(ones(length(male)),male)
Y = data[:,5]
Z = data[:,6:end]
Zmatlab = cat(Z[:,1],Z[:,2],Z[:,3];dims=3)

θ  = [1.255797, -.4242443,  1.095317, -1.063336,  .0436522];

@views @inline function asclogit(bstart::Vector,Y::Array,X::Array,Z::Array,J::Int64,W::Array=ones(length(Y)))

    ## error checking
    @assert ((!isempty(X) || !isempty(Z)) && !isempty(Y))    "You must supply data to the model"
    @assert (ndims(Y)==1 && size(Y,2)==1)                    "Y must be a 1-D Array"
    @assert (minimum(Y)==1 && maximum(Y)==J) "Y should contain integers numbered consecutively from 1 through J"
    if !isempty(X)
        @assert ndims(X)==2          "X must be a 2-dimensional matrix"
        @assert size(X,1)==size(Y,1) "The 1st dimension of X should equal the number of observations in Y"
    end
    if !isempty(Z)
        @assert ndims(Z)==3          "Z must be a 3-dimensional tensor"
        @assert size(Z,1)==size(Y,1) "The 1st dimension of Z should equal the number of observations in Y"
    end

    K1 = size(X,2)
    K2 = size(Z,2)

    function f(b)
        T = promote_type(promote_type(promote_type(eltype(X),eltype(b)),eltype(Z)),eltype(W))
        num   = zeros(T,size(Y))
        dem   = zeros(T,size(Y))
        temp  = zeros(T,size(Y))
        numer =  ones(T,size(Y,1),J)
        P     = zeros(T,size(Y,1),J)
        ℓ     =  zero(T)
        b2 = b[K1*(J-1)+1:K1*(J-1)+K2]

        for j=1:(J-1)
            temp       .= X*b[(j-1)*K1+1:j*K1] .+ (Z[:,:,j].-Z[:,:,J])*b2
            num        .= (Y.==j).*temp.+num
            dem        .+= exp.(temp)
            numer[:,j] .=  exp.(temp)
        end
        dem.+=1
        P   .=numer./(1 .+ sum(numer;dims=2))

        ℓ = -W'*(num.-log.(dem))
    end

    function g!(G,b)
        T     = promote_type(promote_type(promote_type(eltype(X),eltype(b)),eltype(Z)),eltype(W))
        numer = zeros(T,size(Y,1),J)
        P     = zeros(T,size(Y,1),J)
        numg  = zeros(T,K2)
        demg  = zeros(T,K2)
        b2    = b[K1*(J-1)+1:K1*(J-1)+K2]

        G = Array{T}(undef, length(b))
        for j=1:(J-1)
            numer[:,j] .= exp.( X*b[(j-1)*K1+1:j*K1] .+ (Z[:,:,j].-Z[:,:,J])*b2 )
        end
        P   .=numer./(1 .+ sum(numer;dims=2))

        for j=1:(J-1)
            G[(j-1)*K1+1:j*K1] .= -X'*(W.*((Y.==j).-P[:,j]))
        end

        for j=1:(J-1)
            numg .+= -(Z[:,:,j]-Z[:,:,J])'*(W.*(Y.==j))
            demg .+= -(Z[:,:,j]-Z[:,:,J])'*(W.*P[:,j])
        end
        G[K1*(J-1)+1:K1*(J-1)+K2] .= numg.-demg
        G
    end

    td = TwiceDifferentiable(f, bstart, autodiff = :forwarddiff)
    rs = optimize(td, bstart, LBFGS(; linesearch = LineSearches.BackTracking()), Optim.Options(iterations=100_000,g_tol=1e-8,f_tol=1e-8))
    β  = Optim.minimizer(rs)
    ℓ  = Optim.minimum(rs)*(-1)
    H  = Optim.hessian!(td, β)
    se = sqrt.(diag(inv(H)))
    q = Array{Float64}(undef, length(β))
    grad = g!(q,β)
    gradzero = g!(q,zeros(size(β)))

    return β,se,ℓ,grad,gradzero
end

β,se_β,like,grad = asclogit(zeros(size(θ)),Y,X,Zmatlab,3)
β,se_β,like,g1,g0 = @time asclogit(β,Y,X,Zmatlab,3)
println([g1 g0])

## test_autochoice_mew2.jl
using DelimitedFiles
using Base.Iterators: flatten
using LinearAlgebra: diag
using Optim
using DataFrames
using HTTP
using uCSV
using LineSearches

# Read in auto data from Stata example
fname = "https://gist.githubusercontent.com/tyleransom/696f3ef806fa4c84e1902f34b8955ca4/raw/49d3f51abe9f126e48907522a69e14b2fd7b5498/autochoice.csv"
data = DataFrame(uCSV.read(IOBuffer(HTTP.get(fname).body), quotes='"', header=1))
data = Matrix(data)
id = data[:,1]
male = data[:,2]
income = data[:,3]
hhsize = data[:,4]
X = hcat(ones(length(male)),male)
Y = data[:,5]
Z = data[:,6:end]
Zmatlab = cat(Z[:,1],Z[:,2],Z[:,3];dims=3)

θ  = [1.255797, -.4242443,  1.095317, -1.063336,  .0436522];

@views @inline function asclogit(bstart::Vector,Y::Array,X::Array,Z::Array,J::Int64,W::Array=ones(length(Y)))

    ## error checking
    @assert ((!isempty(X) || !isempty(Z)) && !isempty(Y))    "You must supply data to the model"
    @assert (ndims(Y)==1 && size(Y,2)==1)                    "Y must be a 1-D Array"
    @assert (minimum(Y)==1 && maximum(Y)==J) "Y should contain integers numbered consecutively from 1 through J"
    if !isempty(X)
        @assert ndims(X)==2          "X must be a 2-dimensional matrix"
        @assert size(X,1)==size(Y,1) "The 1st dimension of X should equal the number of observations in Y"
    end
    if !isempty(Z)
        @assert ndims(Z)==3          "Z must be a 3-dimensional tensor"
        @assert size(Z,1)==size(Y,1) "The 1st dimension of Z should equal the number of observations in Y"
    end

    K1 = size(X,2)
    K2 = size(Z,2)

    function f(b)
        T = promote_type(promote_type(promote_type(eltype(X),eltype(b)),eltype(Z)),eltype(W))
        num   = zeros(T,size(Y))
        dem   = zeros(T,size(Y))
        temp  = zeros(T,size(Y))
        numer =  ones(T,size(Y,1),J)
        P     = zeros(T,size(Y,1),J)
        ℓ     =  zero(T)
        b2 = b[K1*(J-1)+1:K1*(J-1)+K2]

        for j=1:(J-1)
            temp       .= X*b[(j-1)*K1+1:j*K1] .+ (Z[:,:,j].-Z[:,:,J])*b2
            num        .= (Y.==j).*temp.+num
            dem        .+= exp.(temp)
            numer[:,j] .=  exp.(temp)
        end
        dem.+=1
        P   .=numer./(1 .+ sum(numer;dims=2))

        ℓ = -W'*(num.-log.(dem))
    end

    function g!(G,b)
        T     = promote_type(promote_type(promote_type(eltype(X),eltype(b)),eltype(Z)),eltype(W))
        numer = zeros(T,size(Y,1),J)
        P     = zeros(T,size(Y,1),J)
        numg  = zeros(T,K2)
        demg  = zeros(T,K2)
        b2    = b[K1*(J-1)+1:K1*(J-1)+K2]

        #G = Array{T}(undef, length(b))
        G .= T(0)
        for j=1:(J-1)
            numer[:,j] .= exp.( X*b[(j-1)*K1+1:j*K1] .+ (Z[:,:,j].-Z[:,:,J])*b2 )
        end
        P   .=numer./(1 .+ sum(numer;dims=2))

        for j=1:(J-1)
            G[(j-1)*K1+1:j*K1] .= -X'*(W.*((Y.==j).-P[:,j]))
        end

        for j=1:(J-1)
            numg .+= -(Z[:,:,j]-Z[:,:,J])'*(W.*(Y.==j))
            demg .+= -(Z[:,:,j]-Z[:,:,J])'*(W.*P[:,j])
        end
        G[K1*(J-1)+1:K1*(J-1)+K2] .= numg.-demg
        G
    end

    td = TwiceDifferentiable(f, g!, bstart, autodiff = :forwarddiff)
    rs = optimize(td, bstart, LBFGS(; linesearch = LineSearches.BackTracking()), Optim.Options(iterations=100_000,g_tol=1e-8,f_tol=1e-8))
    β  = Optim.minimizer(rs)
    ℓ  = Optim.minimum(rs)*(-1)
    H  = Optim.hessian!(td, β)
    se = sqrt.(diag(inv(H)))
    q = Array{Float64}(undef, length(β))
    grad = g!(q,β)
    gradzero = g!(q,zeros(size(β)))

    return β,se,ℓ,grad,gradzero
end

β,se_β,like,grad = asclogit(zeros(size(θ)),Y,X,Zmatlab,3)
β,se_β,like,g1,g0 = @time asclogit(β,Y,X,Zmatlab,3)
println([g1 g0])

## test_autochoice_mew3.jl
using DelimitedFiles
using Base.Iterators: flatten
using LinearAlgebra: diag
using Optim
using DataFrames
using HTTP
using uCSV
using LineSearches
using ForwardDiff

# Read in auto data from Stata example
fname = "https://gist.githubusercontent.com/tyleransom/696f3ef806fa4c84e1902f34b8955ca4/raw/49d3f51abe9f126e48907522a69e14b2fd7b5498/autochoice.csv"
data = DataFrame(uCSV.read(IOBuffer(HTTP.get(fname).body), quotes='"', header=1))
data = Matrix(data)
id = data[:,1]
male = data[:,2]
income = data[:,3]
hhsize = data[:,4]
X = hcat(ones(length(male)),male)
Y = data[:,5]
Z = data[:,6:end]
Zmatlab = cat(Z[:,1],Z[:,2],Z[:,3];dims=3)

θ  = [1.255797, -.4242443,  1.095317, -1.063336,  .0436522];

@views @inline function asclogit2(bstart::Vector,Y::Array,X::Array,Z::Array,J::Int64,W::Array=ones(length(Y)))

    ## error checking
    @assert ((!isempty(X) || !isempty(Z)) && !isempty(Y))    "You must supply data to the model"
    @assert (ndims(Y)==1 && size(Y,2)==1)                    "Y must be a 1-D Array"
    @assert (minimum(Y)==1 && maximum(Y)==J) "Y should contain integers numbered consecutively from 1 through J"
    if !isempty(X)
        @assert ndims(X)==2          "X must be a 2-dimensional matrix"
        @assert size(X,1)==size(Y,1) "The 1st dimension of X should equal the number of observations in Y"
    end
    if !isempty(Z)
        @assert ndims(Z)==3          "Z must be a 3-dimensional tensor"
        @assert size(Z,1)==size(Y,1) "The 1st dimension of Z should equal the number of observations in Y"
    end

    K1 = size(X,2)
    K2 = size(Z,2)

    function fg!(F,G,b)
        T = promote_type(promote_type(promote_type(eltype(X),eltype(b)),eltype(Z)),eltype(W))
        num   = zeros(T,size(Y))
        dem   = zeros(T,size(Y))
        temp  = zeros(T,size(Y))
        numer =  ones(T,size(Y,1),J)
        P     = zeros(T,size(Y,1),J)
        b2 = b[K1*(J-1)+1:K1*(J-1)+K2]

        for j=1:(J-1)
            temp       .= X*b[(j-1)*K1+1:j*K1] .+ (Z[:,:,j].-Z[:,:,J])*b2
            num        .= (Y.==j).*temp.+num
            dem        .+= exp.(temp)
            numer[:,j] .=  exp.(temp)
        end
        dem.+=1
        P   .=numer./(1 .+ sum(numer;dims=2))

        if G != nothing
            G .= T(0)
            numg  = zeros(T,K2)
            demg  = zeros(T,K2)
            for j=1:(J-1)
                G[(j-1)*K1+1:j*K1] .= -X'*(W.*((Y.==j).-P[:,j]))
            end
            for j=1:(J-1)
                numg .-= (Z[:,:,j]-Z[:,:,J])'*(W.*(Y.==j))
                demg .-= (Z[:,:,j]-Z[:,:,J])'*(W.*P[:,j])
            end
            G[K1*(J-1)+1:K1*(J-1)+K2] .= numg.-demg
            return nothing
        end

        if F != nothing
            ℓ = T(0)
            ℓ = -W'*(num.-log.(dem))
        end
    end

    od = OnceDifferentiable(Optim.only_fg!(fg!), bstart)
    td = TwiceDifferentiable(od, autodiff = :forwarddiff)
    rs = optimize(td, bstart, LBFGS(; linesearch = LineSearches.BackTracking()), Optim.Options(iterations=100_000,g_tol=1e-6,f_tol=1e-6))
    β  = Optim.minimizer(rs)
    ℓ  = Optim.minimum(rs)*(-1)
    H  = Optim.hessian!(td, β)
    se = sqrt.(diag(inv(H)))
    gf = Array{Float64}(undef, length(β))
    fg!(nothing,gf,β)

    return β,se,ℓ,gf
end

β,se_β,like,grad = @time asclogit2(2*ones(size(θ)),Y,X,Zmatlab,3)
id	male	income	size	choicecat	dealer1	dealer2	dealer3
1	1	46.7	2	3	18	8	5
2	1	26.1	3	1	17	6	2
3	1	32.7	4	1	12	6	2
4	0	49.2	2	2	18	7	4
5	1	24.3	4	1	10	7	3
6	0	39	4	1	13	4	1
7	1	33	2	1	24	7	4
8	1	20.3	4	1	9	3	3
9	1	38	3	2	14	10	2
10	0	60.4	3	1	18	5	4
11	1	69	3	1	23	8	2
12	1	27.7	1	2	21	9	5
13	1	41	2	2	20	10	5
14	1	65.6	4	2	13	4	5
15	1	24.8	3	1	21	5	3
16	1	21.6	2	1	19	8	3
17	0	43.3	1	1	23	9	5
18	0	44.4	1	2	22	10	4
19	1	45.3	4	1	17	8	1
20	1	48.8	3	1	20	10	5
21	1	44.7	4	3	23	7	2
22	1	57.5	1	1	24	9	4
23	1	42	2	1	20	8	3
24	1	40.5	1	1	23	10	4
25	1	44.4	1	1	24	10	5
26	0	50.3	2	1	21	7	4
27	0	41.1	4	1	13	3	2
28	1	49	3	2	22	10	3
29	0	43	3	1	21	7	2
30	1	43.8	4	3	16	8	3
31	1	46.6	2	3	24	8	5
32	1	45.7	3	1	14	10	2
33	1	69.8	3	1	23	10	5
34	0	45.6	2	2	21	7	5
35	1	20.9	1	1	23	10	4
36	1	30.6	4	1	20	4	4
37	0	40.7	2	3	18	7	3
38	1	46.7	2	3	18	8	5
39	1	26.1	3	1	17	6	2
40	1	32.7	4	1	12	6	2
41	0	49.2	2	2	18	7	4
42	1	24.3	4	1	10	7	3
43	0	39	4	1	13	4	1
44	1	33	2	1	24	7	4
45	1	20.3	4	1	9	3	3
46	1	38	3	2	14	10	2
47	0	60.4	3	1	18	5	4
48	1	69	3	1	23	8	2
49	1	27.7	1	2	21	9	5
50	1	41	2	2	20	10	5
51	1	65.6	4	2	13	4	5
52	1	24.8	3	1	21	5	3
53	1	21.6	2	1	19	8	3
54	0	43.3	1	1	23	9	5
55	0	44.4	1	2	22	10	4
56	1	45.3	4	1	17	8	1
57	1	48.8	3	1	20	10	5
58	1	44.7	4	3	23	7	2
59	1	57.5	1	1	24	9	4
60	1	42	2	1	20	8	3
61	1	40.5	1	1	23	10	4
62	1	44.4	1	1	24	10	5
63	0	50.3	2	1	21	7	4
64	0	41.1	4	1	13	3	2
65	1	49	3	2	22	10	3
66	0	43	3	1	21	7	2
67	1	43.8	4	3	16	8	3
68	1	46.6	2	3	24	8	5
69	1	45.7	3	1	14	10	2
70	1	69.8	3	1	23	10	5
71	0	45.6	2	2	21	7	5
72	1	20.9	1	1	23	10	4
73	1	30.6	4	1	20	4	4
74	0	40.7	2	3	18	7	3
75	1	46.7	2	3	18	8	5
76	1	26.1	3	1	17	6	2
77	1	32.7	4	1	12	6	2
78	0	49.2	2	2	18	7	4
79	1	24.3	4	1	10	7	3
80	0	39	4	1	13	4	1
81	1	33	2	1	24	7	4
82	1	20.3	4	1	9	3	3
83	1	38	3	2	14	10	2
84	0	60.4	3	1	18	5	4
85	1	69	3	1	23	8	2
86	1	27.7	1	2	21	9	5
87	1	41	2	2	20	10	5
88	1	65.6	4	2	13	4	5
89	1	24.8	3	1	21	5	3
90	1	21.6	2	1	19	8	3
91	0	43.3	1	1	23	9	5
92	0	44.4	1	2	22	10	4
93	1	45.3	4	1	17	8	1
94	1	48.8	3	1	20	10	5
95	1	44.7	4	3	23	7	2
96	1	57.5	1	1	24	9	4
97	1	42	2	1	20	8	3
98	1	40.5	1	1	23	10	4
99	1	44.4	1	1	24	10	5
100	0	50.3	2	1	21	7	4
101	0	41.1	4	1	13	3	2
102	1	49	3	2	22	10	3
103	0	43	3	1	21	7	2
104	1	43.8	4	3	16	8	3
105	1	46.6	2	3	24	8	5
106	1	45.7	3	1	14	10	2
107	1	69.8	3	1	23	10	5
108	0	45.6	2	2	21	7	5
109	1	20.9	1	1	23	10	4
110	1	30.6	4	1	20	4	4
111	0	40.7	2	3	18	7	3
112	1	46.7	2	3	18	8	5
113	1	26.1	3	1	17	6	2
114	1	32.7	4	1	12	6	2
115	0	49.2	2	2	18	7	4
116	1	24.3	4	1	10	7	3
117	0	39	4	1	13	4	1
118	1	33	2	1	24	7	4
119	1	20.3	4	1	9	3	3
120	1	38	3	2	14	10	2
121	0	60.4	3	1	18	5	4
122	1	69	3	1	23	8	2
123	1	27.7	1	2	21	9	5
124	1	41	2	2	20	10	5
125	1	65.6	4	2	13	4	5
126	1	24.8	3	1	21	5	3
127	1	21.6	2	1	19	8	3
128	0	43.3	1	1	23	9	5
129	0	44.4	1	2	22	10	4
130	1	45.3	4	1	17	8	1
131	1	48.8	3	1	20	10	5
132	1	44.7	4	3	23	7	2
133	1	57.5	1	1	24	9	4
134	1	42	2	1	20	8	3
135	1	40.5	1	1	23	10	4
136	1	44.4	1	1	24	10	5
137	0	50.3	2	1	21	7	4
138	0	41.1	4	1	13	3	2
139	1	49	3	2	22	10	3
140	0	43	3	1	21	7	2
141	1	43.8	4	3	16	8	3
142	1	46.6	2	3	24	8	5
143	1	45.7	3	1	14	10	2
144	1	69.8	3	1	23	10	5
145	0	45.6	2	2	21	7	5
146	1	20.9	1	1	23	10	4
147	1	30.6	4	1	20	4	4
148	0	40.7	2	3	18	7	3
149	1	46.7	2	3	18	8	5
150	1	26.1	3	1	17	6	2
151	1	32.7	4	1	12	6	2
152	0	49.2	2	2	18	7	4
153	1	24.3	4	1	10	7	3
154	0	39	4	1	13	4	1
155	1	33	2	1	24	7	4
156	1	20.3	4	1	9	3	3
157	1	38	3	2	14	10	2
158	0	60.4	3	1	18	5	4
159	1	69	3	1	23	8	2
160	1	27.7	1	2	21	9	5
161	1	41	2	2	20	10	5
162	1	65.6	4	2	13	4	5
163	1	24.8	3	1	21	5	3
164	1	21.6	2	1	19	8	3
165	0	43.3	1	1	23	9	5
166	0	44.4	1	2	22	10	4
167	1	45.3	4	1	17	8	1
168	1	48.8	3	1	20	10	5
169	1	44.7	4	3	23	7	2
170	1	57.5	1	1	24	9	4
171	1	42	2	1	20	8	3
172	1	40.5	1	1	23	10	4
173	1	44.4	1	1	24	10	5
174	0	50.3	2	1	21	7	4
175	0	41.1	4	1	13	3	2
176	1	49	3	2	22	10	3
177	0	43	3	1	21	7	2
178	1	43.8	4	3	16	8	3
179	1	46.6	2	3	24	8	5
180	1	45.7	3	1	14	10	2
181	1	69.8	3	1	23	10	5
182	0	45.6	2	2	21	7	5
183	1	20.9	1	1	23	10	4
184	1	30.6	4	1	20	4	4
185	0	40.7	2	3	18	7	3
186	1	46.7	2	3	18	8	5
187	1	26.1	3	1	17	6	2
188	1	32.7	4	1	12	6	2
189	0	49.2	2	2	18	7	4
190	1	24.3	4	1	10	7	3
191	0	39	4	1	13	4	1
192	1	33	2	1	24	7	4
193	1	20.3	4	1	9	3	3
194	1	38	3	2	14	10	2
195	0	60.4	3	1	18	5	4
196	1	69	3	1	23	8	2
197	1	27.7	1	2	21	9	5
198	1	41	2	2	20	10	5
199	1	65.6	4	2	13	4	5
200	1	24.8	3	1	21	5	3
201	1	21.6	2	1	19	8	3
202	0	43.3	1	1	23	9	5
203	0	44.4	1	2	22	10	4
204	1	45.3	4	1	17	8	1
205	1	48.8	3	1	20	10	5
206	1	44.7	4	3	23	7	2
207	1	57.5	1	1	24	9	4
208	1	42	2	1	20	8	3
209	1	40.5	1	1	23	10	4
210	1	44.4	1	1	24	10	5
211	0	50.3	2	1	21	7	4
212	0	41.1	4	1	13	3	2
213	1	49	3	2	22	10	3
214	0	43	3	1	21	7	2
215	1	43.8	4	3	16	8	3
216	1	46.6	2	3	24	8	5
217	1	45.7	3	1	14	10	2
218	1	69.8	3	1	23	10	5
219	0	45.6	2	2	21	7	5
220	1	20.9	1	1	23	10	4
221	1	30.6	4	1	20	4	4
222	0	40.7	2	3	18	7	3
223	1	46.7	2	3	18	8	5
224	1	26.1	3	1	17	6	2
225	1	32.7	4	1	12	6	2
226	0	49.2	2	2	18	7	4
227	1	24.3	4	1	10	7	3
228	0	39	4	1	13	4	1
229	1	33	2	1	24	7	4
230	1	20.3	4	1	9	3	3
231	1	38	3	2	14	10	2
232	0	60.4	3	1	18	5	4
233	1	69	3	1	23	8	2
234	1	27.7	1	2	21	9	5
235	1	41	2	2	20	10	5
236	1	65.6	4	2	13	4	5
237	1	24.8	3	1	21	5	3
238	1	21.6	2	1	19	8	3
239	0	43.3	1	1	23	9	5
240	0	44.4	1	2	22	10	4
241	1	45.3	4	1	17	8	1
242	1	48.8	3	1	20	10	5
243	1	44.7	4	3	23	7	2
244	1	57.5	1	1	24	9	4
245	1	42	2	1	20	8	3
246	1	40.5	1	1	23	10	4
247	1	44.4	1	1	24	10	5
248	0	50.3	2	1	21	7	4
249	0	41.1	4	1	13	3	2
250	1	49	3	2	22	10	3
251	0	43	3	1	21	7	2
252	1	43.8	4	3	16	8	3
253	1	46.6	2	3	24	8	5
254	1	45.7	3	1	14	10	2
255	1	69.8	3	1	23	10	5
256	0	45.6	2	2	21	7	5
257	1	20.9	1	1	23	10	4
258	1	30.6	4	1	20	4	4
259	0	40.7	2	3	18	7	3
260	1	46.7	2	3	18	8	5
261	1	26.1	3	1	17	6	2
262	1	32.7	4	1	12	6	2
263	0	49.2	2	2	18	7	4
264	1	24.3	4	1	10	7	3
265	0	39	4	1	13	4	1
266	1	33	2	1	24	7	4
267	1	20.3	4	1	9	3	3
268	1	38	3	2	14	10	2
269	0	60.4	3	1	18	5	4
270	1	69	3	1	23	8	2
271	1	27.7	1	2	21	9	5
272	1	41	2	2	20	10	5
273	1	65.6	4	2	13	4	5
274	1	24.8	3	1	21	5	3
275	1	21.6	2	1	19	8	3
276	0	43.3	1	1	23	9	5
277	0	44.4	1	2	22	10	4
278	1	45.3	4	1	17	8	1
279	1	48.8	3	1	20	10	5
280	1	44.7	4	3	23	7	2
281	1	57.5	1	1	24	9	4
282	1	42	2	1	20	8	3
283	1	40.5	1	1	23	10	4
284	1	44.4	1	1	24	10	5
285	0	50.3	2	1	21	7	4
286	0	41.1	4	1	13	3	2
287	1	49	3	2	22	10	3
288	0	43	3	1	21	7	2
289	1	43.8	4	3	16	8	3
290	1	46.6	2	3	24	8	5
291	1	45.7	3	1	14	10	2
292	1	69.8	3	1	23	10	5
293	0	45.6	2	2	21	7	5
294	1	20.9	1	1	23	10	4
295	1	30.6	4	1	20	4	4
	using DelimitedFiles
	using Base.Iterators: flatten
	using LinearAlgebra: diag
	using Optim
	using DataFrames
	using HTTP
	using uCSV
	using LineSearches

	# Read in auto data from Stata example
	fname = "https://gist.githubusercontent.com/tyleransom/696f3ef806fa4c84e1902f34b8955ca4/raw/49d3f51abe9f126e48907522a69e14b2fd7b5498/autochoice.csv"
	data = DataFrame(uCSV.read(IOBuffer(HTTP.get(fname).body), quotes='"', header=1))
	data = Matrix(data)
	id = data[:,1]
	male = data[:,2]
	income = data[:,3]
	hhsize = data[:,4]
	X = hcat(ones(length(male)),male)
	Y = data[:,5]
	Z = data[:,6:end]
	Zmatlab = cat(Z[:,1],Z[:,2],Z[:,3];dims=3)

	θ = [1.255797, -.4242443, 1.095317, -1.063336, .0436522];

	@views @inline function asclogit(bstart::Vector,Y::Array,X::Array,Z::Array,J::Int64,W::Array=ones(length(Y)))

	## error checking
	@assert ((!isempty(X) \|\| !isempty(Z)) && !isempty(Y)) "You must supply data to the model"
	@assert (ndims(Y)==1 && size(Y,2)==1) "Y must be a 1-D Array"
	@assert (minimum(Y)==1 && maximum(Y)==J) "Y should contain integers numbered consecutively from 1 through J"
	if !isempty(X)
	@assert ndims(X)==2 "X must be a 2-dimensional matrix"
	@assert size(X,1)==size(Y,1) "The 1st dimension of X should equal the number of observations in Y"
	end
	if !isempty(Z)
	@assert ndims(Z)==3 "Z must be a 3-dimensional tensor"
	@assert size(Z,1)==size(Y,1) "The 1st dimension of Z should equal the number of observations in Y"
	end

	K1 = size(X,2)
	K2 = size(Z,2)

	function f(b)
	T = promote_type(promote_type(promote_type(eltype(X),eltype(b)),eltype(Z)),eltype(W))
	num = zeros(T,size(Y))
	dem = zeros(T,size(Y))
	temp = zeros(T,size(Y))
	numer = ones(T,size(Y,1),J)
	P = zeros(T,size(Y,1),J)
	ℓ = zero(T)
	b2 = b[K1(J-1)+1:K1(J-1)+K2]

	for j=1:(J-1)
	temp .= Xb[(j-1)K1+1:jK1] .+ (Z[:,:,j].-Z[:,:,J])b2
	num .= (Y.==j).*temp.+num
	dem .+= exp.(temp)
	numer[:,j] .= exp.(temp)
	end
	dem.+=1
	P .=numer./(1 .+ sum(numer;dims=2))

	ℓ = -W'*(num.-log.(dem))
	end

	function g!(G,b)
	T = promote_type(promote_type(promote_type(eltype(X),eltype(b)),eltype(Z)),eltype(W))
	numer = zeros(T,size(Y,1),J)
	P = zeros(T,size(Y,1),J)
	numg = zeros(T,K2)
	demg = zeros(T,K2)
	b2 = b[K1(J-1)+1:K1(J-1)+K2]

	G = Array{T}(undef, length(b))
	for j=1:(J-1)
	numer[:,j] .= exp.( Xb[(j-1)K1+1:jK1] .+ (Z[:,:,j].-Z[:,:,J])b2 )
	end
	P .=numer./(1 .+ sum(numer;dims=2))

	for j=1:(J-1)
	G[(j-1)K1+1:jK1] .= -X'(W.((Y.==j).-P[:,j]))
	end

	for j=1:(J-1)
	numg .+= -(Z[:,:,j]-Z[:,:,J])'(W.(Y.==j))
	demg .+= -(Z[:,:,j]-Z[:,:,J])'(W.P[:,j])
	end
	G[K1(J-1)+1:K1(J-1)+K2] .= numg.-demg
	G
	end

	td = TwiceDifferentiable(f, bstart, autodiff = :forwarddiff)
	rs = optimize(td, bstart, LBFGS(; linesearch = LineSearches.BackTracking()), Optim.Options(iterations=100_000,g_tol=1e-8,f_tol=1e-8))
	β = Optim.minimizer(rs)
	ℓ = Optim.minimum(rs)*(-1)
	H = Optim.hessian!(td, β)
	se = sqrt.(diag(inv(H)))
	q = Array{Float64}(undef, length(β))
	grad = g!(q,β)
	gradzero = g!(q,zeros(size(β)))

	return β,se,ℓ,grad,gradzero
	end

	β,se_β,like,grad = asclogit(zeros(size(θ)),Y,X,Zmatlab,3)
	β,se_β,like,g1,g0 = @time asclogit(β,Y,X,Zmatlab,3)
	println([g1 g0])