text:
function Y_l = load_regressor_UR10E(in1,in2,in3)
%LOAD_REGRESSOR_UR10E
%    Y_L = LOAD_REGRESSOR_UR10E(IN1,IN2,IN3)

%    This function was generated by the Symbolic Math Toolbox version 8.2.
%    20-Oct-2021 07:53:37

q2 = in1(2,:);
q3 = in1(3,:);
q4 = in1(4,:);
q5 = in1(5,:);
q6 = in1(6,:);
q2d1 = in3(1,:);
q2d2 = in3(2,:);
q2d3 = in3(3,:);
q2d4 = in3(4,:);
q2d5 = in3(5,:);
q2d6 = in3(6,:);
qd1 = in2(1,:);
qd2 = in2(2,:);
qd3 = in2(3,:);
qd4 = in2(4,:);
qd5 = in2(5,:);
qd6 = in2(6,:);
t2 = sin(1.57079632679);
t3 = cos(q2);
t4 = sin(q3);
t5 = cos(q3);
t6 = sin(q2);
t7 = sin(q5);
t8 = cos(q4);
t9 = sin(q4);
t10 = t2.^2;
t11 = qd1.*t2.*t3.*t5;
t24 = qd1.*t2.*t4.*t6;
t12 = t11-t24;
t13 = qd1.*t2.*t3.*t4;
t14 = qd1.*t2.*t5.*t6;
t15 = t13+t14;
t16 = t2.*t3.*t4;
t17 = t2.*t5.*t6;
t18 = t16+t17;
t19 = t2.*t3.*t5;
t21 = t2.*t4.*t6;
t20 = t19-t21;
t22 = t2.*t9.*t18;
t23 = t2.*t8.*t15;
t25 = t2.*t9.*t12;
t26 = t23+t25;
t27 = cos(q5);
t28 = qd2+qd3+qd4;
t29 = t27.*t28;
t41 = t7.*t26;
t30 = qd6+t29-t41;
t40 = t2.*t8.*t20;
t31 = t22-t40;
t32 = t7.*t10.*t30.*t31;
t33 = t7.^2;
t34 = t2.*t8.*t12;
t42 = t2.*t9.*t15;
t35 = t34-t42;
t36 = t2.*t8.*t18;
t37 = t2.*t9.*t20;
t38 = t36+t37;
t39 = t10.*t33.*t35.*t38;
t43 = sin(q6);
t44 = cos(q6);
t45 = t26.*t27;
t46 = t7.*t28;
t47 = t45+t46;
t48 = t38.*t43;
t49 = t26.*t43;
t62 = t27.*t35.*t44;
t50 = t49-t62;
t51 = t7.*t10.*t38.*t50;
t52 = qd5-t34+t42;
t53 = t44.*t47;
t63 = t43.*t52;
t54 = t53-t63;
t55 = t7.*t10.*t31.*t54;
t56 = t31.*t43;
t64 = t27.*t38.*t44;
t57 = t56-t64;
t58 = t7.*t10.*t35.*t57;
t59 = t27.*t31.*t44;
t60 = t48+t59;
t65 = t10.*t30.*t60;
t61 = t51+t55+t58-t65;
t66 = t10.*t27.*t57;
t67 = t10.*t33.*t38.*t44;
t68 = t31.*t44;
t69 = t27.*t38.*t43;
t70 = t68+t69;
t71 = t2.*t27.*t70;
t72 = t44.*t52;
t73 = t43.*t47;
t74 = t72+t73;
t75 = t29-t41;
t76 = t38.*t44;
t84 = t27.*t31.*t43;
t77 = t76-t84;
t78 = t2.*t30.*t77;
t79 = t26.*t44;
t80 = t27.*t35.*t43;
t81 = t79+t80;
t82 = t2.*t7.*t74.*(t22-t40);
t85 = t2.*t7.*t38.*t81;
t86 = t2.*t7.*t35.*t70;
t83 = t78+t82-t85-t86;
t87 = t10.*t50.*t57;
t88 = t27.*t44.*(t22-t40);
t89 = t48+t88;
t90 = t56-t64;
t91 = t2.*t7.*t43.*t57;
t93 = t2.*t7.*t44.*t70;
t92 = t91-t93;
t94 = t2.*t57.*t81;
t95 = t2.*t50.*t70;
t96 = t2.*t74.*t89;
t97 = t74.*t77;
t98 = t81.*(t68+t69);
t99 = t97+t98;
t100 = t68+t69;
t101 = t2.*t3.*t5.*(3.9e1./1.0e3);
t109 = t2.*t4.*t6.*(3.9e1./1.0e3);
t102 = t101-t109;
t103 = t2.*t3.*t4.*(3.9e1./1.0e3);
t104 = t2.*t5.*t6.*(3.9e1./1.0e3);
t105 = t103+t104;
t106 = t2.*t8.*t18.*(2.7e1./2.0e2);
t107 = t2.*t8.*t105;
t108 = t2.*t9.*t20.*(2.7e1./2.0e2);
t110 = t2.*t9.*t102;
t111 = t106+t107+t108+t110;
t112 = t2.*t9.*t18.*(2.7e1./2.0e2);
t113 = t2.*t9.*t105;
t120 = t2.*t8.*t20.*(2.7e1./2.0e2);
t121 = t2.*t8.*t102;
t114 = t112+t113-t120-t121;
t115 = t27.*t114;
t116 = t2.*t3.*(6.13e2./1.0e3);
t117 = t2.*t3.*t5.*(5.71e2./1.0e3);
t122 = t2.*t4.*t6.*(5.71e2./1.0e3);
t118 = t116+t117-t122;
t119 = t7.*t38.*(3.0./2.5e1);
t123 = t7.*t114;
t124 = t27.*t118;
t244 = t27.*t38.*(3.0./2.5e1);
t125 = t123+t124-t244;
t126 = qd2.*(6.13e2./1.0e3);
t127 = qd1.*t2.*t6.*(3.9e1./1.0e3);
t128 = t126+t127;
t129 = qd1.*t2.*t3.*(6.13e2./1.0e3);
t130 = qd1.*t2.*t3.*t5.*(5.71e2./1.0e3);
t155 = qd1.*t2.*t4.*t6.*(5.71e2./1.0e3);
t131 = t129+t130-t155;
t132 = t4.*t128;
t154 = qd1.*t2.*t3.*t5.*(3.9e1./1.0e3);
t133 = t132-t154;
t134 = t2.*t8.*t133;
t135 = t2.*t9.*t15.*(2.7e1./2.0e2);
t136 = qd2.*(5.71e2./1.0e3);
t137 = qd3.*(5.71e2./1.0e3);
t138 = t5.*t128;
t139 = qd1.*t2.*t3.*t4.*(3.9e1./1.0e3);
t140 = t136+t137+t138+t139;
t141 = t2.*t9.*t140;
t156 = t2.*t8.*t12.*(2.7e1./2.0e2);
t142 = t134+t135+t141-t156;
t143 = t26.*t27.*(3.0./2.5e1);
t144 = t7.*t28.*(3.0./2.5e1);
t149 = t7.*t118;
t145 = t115+t119-t149;
t146 = t43.*t145;
t147 = t31.*t43.*(1.17e2./1.0e3);
t148 = t43.*t111;
t150 = t44.*(t22-t40).*(1.17e2./1.0e3);
t151 = t27.*t38.*t43.*(1.17e2./1.0e3);
t152 = t44.*t145;
t153 = t148+t150+t151+t152;
t157 = t2.*t8.*t15.*(2.7e1./2.0e2);
t158 = t2.*t9.*t12.*(2.7e1./2.0e2);
t159 = t2.*t8.*t140;
t192 = t2.*t9.*t133;
t160 = t157+t158+t159-t192;
t161 = t27.*t142;
t162 = t7.*t26.*(3.0./2.5e1);
t249 = t27.*t131;
t250 = t7.*t142;
t163 = t143+t144-t249-t250;
t201 = t44.*t111;
t202 = t27.*t38.*t44.*(1.17e2./1.0e3);
t164 = t146+t147-t201-t202;
t194 = t7.*t131;
t195 = t27.*t28.*(3.0./2.5e1);
t165 = t161+t162-t194-t195;
t166 = t27.*t111;
t167 = t2.*t6.*(6.13e2./1.0e3);
t168 = t2.*t3.*t4.*(5.71e2./1.0e3);
t169 = t2.*t5.*t6.*(5.71e2./1.0e3);
t170 = t167+t168+t169;
t171 = t7.*t170;
t215 = t7.*t31.*(3.0./2.5e1);
t172 = t166+t171-t215;
t183 = qd1.*t2.*t4.*t6.*(3.9e1./1.0e3);
t173 = t154-t183;
t174 = qd1.*t2.*t5.*t6.*(3.9e1./1.0e3);
t175 = t139+t174;
t176 = qd1.*t2.*t6.*(6.13e2./1.0e3);
t177 = qd1.*t2.*t3.*t4.*(5.71e2./1.0e3);
t178 = qd1.*t2.*t5.*t6.*(5.71e2./1.0e3);
t179 = t176+t177+t178;
t180 = t7.*t179;
t181 = t7.*t35.*(3.0./2.5e1);
t182 = t2.*t8.*t175;
t184 = t2.*t9.*t173;
t185 = t157+t158+t182+t184;
t186 = t27.*t185;
t187 = t180+t181+t186;
t188 = t2.*t9.*t175;
t284 = t2.*t8.*t173;
t189 = t135-t156+t188-t284;
t190 = t44.*t52.*(1.17e2./1.0e3);
t191 = t43.*t47.*(1.17e2./1.0e3);
t193 = t43.*t160;
t196 = t44.*t165;
t197 = t190+t191+t193+t196;
t198 = t44.*t47.*(1.17e2./1.0e3);
t199 = t44.*t160;
t222 = t43.*t52.*(1.17e2./1.0e3);
t223 = t43.*t165;
t200 = t198+t199-t222-t223;
t203 = t26.*t43.*(1.17e2./1.0e3);
t204 = t138+t139;
t205 = t177+t178;
t206 = t7.*t205;
t207 = t2.*t8.*t204;
t208 = t157+t158-t192+t207;
t209 = t27.*t208;
t210 = t181+t206+t209;
t211 = t26.*t44.*(1.17e2./1.0e3);
t212 = t2.*t9.*t204;
t213 = t134+t135-t156+t212;
t214 = t27.*t35.*t43.*(1.17e2./1.0e3);
t216 = t38.*t44.*(1.17e2./1.0e3);
t217 = t168+t169;
t218 = t7.*t217;
t219 = t166-t215+t218;
t220 = t38.*t43.*(1.17e2./1.0e3);
t221 = t2.*t77.*t197;
t224 = t2.*t81.*t153;
t225 = t2.*t50.*t164;
t226 = t5.*(6.13e2./1.0e3);
t227 = t226+5.71e2./1.0e3;
t228 = t2.*t4.*t9.*(6.13e2./1.0e3);
t253 = t2.*t8.*t227;
t229 = t228-t253;
t230 = t27.*(3.0./2.5e1);
t231 = t2.*t9.*t227;
t232 = t2.*t4.*t8.*(6.13e2./1.0e3);
t233 = t231+t232;
t254 = t27.*t233;
t234 = t230-t254;
t235 = t44.*t114;
t236 = t166-t215;
t237 = t27.*t31.*t44.*(1.17e2./1.0e3);
t238 = t27.*t160;
t239 = t181+t238;
t240 = t7.*t44.*(1.17e2./1.0e3);
t245 = t2.*t9.*t27.*(5.71e2./1.0e3);
t241 = t230-t245;
t242 = t2.*t7.*t44.*t164;
t243 = t7.*t43.*(1.17e2./1.0e3);
t246 = t43.*t241;
t247 = t2.*t8.*t44.*(5.71e2./1.0e3);
t248 = t240+t246+t247;
t251 = t2.*t27.*t164;
t252 = t7.*(3.0./2.5e1);
t255 = t43.*t234;
t327 = t44.*t229;
t256 = t240+t255-t327;
t257 = t2.*t7.*t43.*t125;
t258 = t43.*t219;
t259 = t44.*(t112+t113-t120-t121);
t260 = t27.*t44.*(t22-t40).*(1.17e2./1.0e3);
t261 = t220+t258+t259+t260;
t262 = t43.*t210;
t263 = t44.*t213;
t278 = t27.*t35.*t44.*(1.17e2./1.0e3);
t264 = t203+t262+t263-t278;
t265 = t27.*t43.*(3.0./2.5e1);
t266 = t240+t265;
t267 = t43.*t125;
t343 = t7.*t38.*t44.*(1.17e2./1.0e3);
t268 = t267-t343;
t269 = t44.*t75.*(1.17e2./1.0e3);
t358 = t43.*t163;
t270 = t269-t358;
t271 = t27.*t35.*(3.0./2.5e1);
t272 = t7.*t111;
t273 = t27.*t31.*(3.0./2.5e1);
t274 = t43.*t236;
t275 = t220+t235+t237+t274;
t276 = t44.*t142;
t277 = t43.*t239;
t279 = t2.*t77.*t163;
t280 = t43.*t172;
t281 = t2.*t7.*t31.*t200;
t282 = t2.*t7.*t35.*t164;
t283 = t43.*t187;
t285 = t44.*t189;
t286 = t27.*t205;
t335 = t7.*t208;
t287 = t271+t286-t335;
t336 = t27.*t217;
t288 = t272+t273-t336;
t289 = t44.*t219;
t311 = t43.*t114;
t313 = t27.*t31.*t43.*(1.17e2./1.0e3);
t290 = t216+t289-t311-t313;
t291 = t44.*t210;
t337 = t43.*t213;
t292 = t211+t214+t291-t337;
t331 = t2.*t7.*t9.*(5.71e2./1.0e3);
t293 = t252-t331;
t294 = t2.*t8.*t43.*(5.71e2./1.0e3);
t332 = t44.*t241;
t295 = t243+t294-t332;
t296 = t44.*t125;
t297 = t7.*t38.*t43.*(1.17e2./1.0e3);
t298 = t296+t297;
t299 = t43.*t75.*(1.17e2./1.0e3);
t300 = t44.*t163;
t301 = t299+t300;
t302 = t10.*t27.*t153;
t328 = t7.*t233;
t303 = t252-t328;
t304 = t7.*t10.*t44.*t125;
t305 = t43.*t229;
t306 = t44.*t234;
t307 = t44.*(t115+t119-t149);
t308 = t148+t150+t151+t307;
t330 = t27.*t44.*(3.0./2.5e1);
t309 = t243-t330;
t310 = t272+t273;
t312 = t44.*t236;
t346 = t7.*t160;
t314 = t271-t346;
t315 = t10.*t89.*t163;
t316 = t44.*t239;
t345 = t43.*t142;
t317 = t211+t214+t316-t345;
t318 = t7.*t10.*t31.*t197;
t319 = t27.*t179;
t340 = t7.*t185;
t320 = t271+t319-t340;
t339 = t27.*t170;
t321 = t272+t273-t339;
t322 = t44.*t172;
t323 = t44.*t187;
t341 = t43.*t189;
t324 = t211+t214+t323-t341;
t325 = t146+t147-t201-t202;
t326 = t123+t124-t244;
t329 = -t243+t305+t306;
t333 = t44.*(t161+t162-t194-t195);
t334 = t190+t191+t193+t333;
t338 = t220+t235+t237+t280;
t342 = t216-t311-t313+t322;
t344 = t203+t276+t277-t278;
t347 = t216-t311+t312-t313;
t348 = t27.^2;
t349 = t10.*t27.*t50;
t350 = t10.*t33.*t35.*t44;
t351 = t349+t350;
t352 = t2.*t27.*t81;
t354 = t2.*t33.*t35.*t43;
t353 = t352-t354;
t355 = t44.^2;
t356 = t2.*t7.*t43.*t50;
t357 = t43.^2;
t359 = t44.*(t135-t156+t188-t284);
t360 = t203-t278+t283+t359;
t361 = t2.*t50.*t256;
t362 = t2.*t81.*t329;
t363 = t2.*t5.*t9.*(6.13e2./1.0e3);
t364 = t232+t363;
t385 = t2.*t5.*t8.*(6.13e2./1.0e3);
t365 = t228-t385;
t366 = t2.*t57.*t256;
t367 = t2.*t7.*t44.*t256;
t368 = t2.*t7.*t38.*t256;
t369 = t251+t257+t368-t2.*t70.*t303;
t370 = t2.*t27.*t256;
t371 = t2.*t7.*t43.*t303;
t372 = t44.*t233;
t412 = t27.*t43.*t229;
t373 = t372-t412;
t374 = t2.*t81.*t303;
t375 = t2.*t3.*t4.*(9.81e2./1.0e2);
t376 = t2.*t5.*t6.*(9.81e2./1.0e2);
t377 = t375+t376;
t378 = t2.*t3.*t5.*(9.81e2./1.0e2);
t380 = t2.*t4.*t6.*(9.81e2./1.0e2);
t379 = t378-t380;
t381 = t2.*t9.*t377;
t382 = t27.*t44.*(1.17e2./1.0e3);
t426 = t43.*t303;
t383 = t382-t426;
t384 = t44.*t364;
t405 = t27.*t43.*t365;
t386 = t384-t405;
t387 = t43.*t364;
t388 = t27.*t43.*(1.17e2./1.0e3);
t389 = t44.*t303;
t390 = t388+t389;
t391 = t10.*t50.*t303;
t392 = t43.*t233;
t393 = t27.*t44.*t229;
t394 = t392+t393;
t395 = t7.*t10.*t35.*t329;
t396 = t10.*t27.*t308;
t397 = t10.*t57.*t303;
t398 = t7.*t10.*t38.*t329;
t399 = t7.*t10.*t44.*t303;
t400 = t2.*t8.*t377;
t401 = t2.*t9.*t379;
t402 = t400+t401;
t421 = t2.*t8.*t379;
t403 = t381-t421;
t404 = -t243+t305+t306;
t406 = t27.*t44.*(t228-t385);
t407 = t387+t406;
t408 = t164.*t256;
t409 = t10.*t125.*t303;
t410 = t10.*t308.*t329;
t411 = t408+t409+t410;
t413 = t2.*t6.*(1.81e2./1.0e3);
t414 = t413-6.13e2./1.0e3;
t415 = t4.*t414;
t422 = t2.*t3.*t5.*(1.81e2./1.0e3);
t416 = t415-t422;
t417 = t5.*t414;
t418 = t2.*t3.*t4.*(1.81e2./1.0e3);
t419 = t417+t418-5.71e2./1.0e3;
t420 = t44.*t402;
t438 = t2.*t4.*t6.*(1.81e2./1.0e3);
t423 = t422-t438;
t424 = t2.*t5.*t6.*(1.81e2./1.0e3);
t425 = t418+t424;
t427 = t7.*(8.7e1./5.0e2);
t428 = t2.*t9.*t416;
t443 = t2.*t8.*t419;
t429 = t428-t443;
t580 = t27.*t429;
t430 = t427-t580;
t431 = t2.*t8.*t416;
t432 = t2.*t9.*t419;
t433 = t431+t432+3.0./2.5e1;
t434 = t43.*t402;
t435 = t27.*t44.*t403;
t436 = t434+t435;
t437 = t2.*t8.*t425;
t439 = t2.*t9.*t423;
t440 = t437+t439;
t441 = t2.*t8.*t423;
t444 = t2.*t9.*t425;
t442 = t441-t444;
t445 = t7.*t10.*t30;
t446 = t10.*t27.*t47;
t447 = t445+t446;
t448 = q2d6.*t10.*t27;
t449 = q2d2.*t10.*t348;
t450 = q2d3.*t10.*t348;
t451 = q2d4.*t10.*t348;
t452 = t7.*t10.*t30.*t35;
t607 = qd5.*t447;
t608 = qd2.*t7.*t10.*t27.*t35;
t609 = qd3.*t7.*t10.*t27.*t35;
t610 = qd4.*t7.*t10.*t27.*t35;
t611 = q2d1.*t7.*t10.*t27.*t38;
t453 = t448+t449+t450+t451+t452-t607-t608-t609-t610-t611;
t454 = qd2.*t351;
t455 = qd3.*t351;
t456 = qd4.*t351;
t457 = t7.*t10.*t54;
t458 = t7.*t10.*t44.*t47;
t612 = t10.*t27.*t30.*t44;
t613 = t10.*t27.*t44.*(t29-t41);
t459 = qd5.*(t457+t458-t612-t613);
t460 = q2d1.*(t66+t67);
t461 = t10.*t27.*t74;
t462 = t7.*t10.*t30.*t43;
t463 = t461+t462;
t464 = qd6.*t463;
t465 = q2d5.*t10.*t27.*t43;
t614 = t10.*t30.*t50;
t615 = t7.*t10.*t35.*t54;
t616 = q2d6.*t7.*t10.*t44;
t617 = q2d2.*t7.*t10.*t27.*t44.*2.0;
t618 = q2d3.*t7.*t10.*t27.*t44.*2.0;
t619 = q2d4.*t7.*t10.*t27.*t44.*2.0;
t466 = t454+t455+t456+t459+t460+t464+t465-t614-t615-t616-t617-t618-t619;
t467 = t2.*t7.*t74;
t468 = t2.*t7.*t43.*t47;
t620 = t2.*t27.*t43.*t75;
t621 = t2.*t27.*t30.*t43;
t469 = t467+t468-t620-t621;
t623 = t2.*t33.*t38.*t43;
t470 = q2d1.*(t71-t623);
t471 = t2.*t27.*t54;
t472 = t2.*t7.*t30.*t44;
t473 = qd6.*(t471+t472);
t474 = qd2.*t353;
t475 = qd3.*t353;
t476 = qd4.*t353;
t477 = t2.*t7.*t35.*t74;
t478 = q2d5.*t2.*t27.*t44;
t479 = q2d6.*t2.*t7.*t43;
t480 = q2d2.*t2.*t7.*t27.*t43.*2.0;
t481 = q2d3.*t2.*t7.*t27.*t43.*2.0;
t482 = q2d4.*t2.*t7.*t27.*t43.*2.0;
t622 = qd5.*t469;
t624 = t2.*t30.*t81;
t483 = t470+t473+t474+t475+t476+t477+t478+t479+t480+t481+t482-t622-t624;
t484 = t7.*t10.*t44.*t74;
t485 = t7.*t10.*t43.*t54;
t486 = t484+t485;
t487 = t7.*t10.*t75.*t355;
t488 = t10.*t27.*t44.*t54;
t489 = qd5.*(t487+t488);
t490 = t10.*t50.*(t53-t63);
t491 = q2d2.*t10.*t33.*t355;
t492 = q2d3.*t10.*t33.*t355;
t493 = q2d4.*t10.*t33.*t355;
t625 = qd6.*t486;
t626 = q2d5.*t7.*t10.*t43.*t44;
t627 = q2d1.*t7.*t10.*t44.*t57;
t628 = qd2.*t7.*t10.*t44.*t50;
t629 = qd3.*t7.*t10.*t44.*t50;
t630 = qd4.*t7.*t10.*t44.*t50;
t494 = t489+t490+t491+t492+t493-t625-t626-t627-t628-t629-t630;
t495 = q2d1.*t92;
t497 = t2.*t7.*t44.*t81;
t496 = qd2.*(t356-t497);
t498 = t2.*t7.*t44.*t54.*2.0;
t633 = t2.*t7.*t43.*t74.*2.0;
t499 = t498-t633;
t500 = t2.*t27.*t44.*t74;
t501 = t2.*t27.*t43.*t54;
t502 = t2.*t7.*t43.*t44.*t75.*2.0;
t503 = t500+t501+t502;
t504 = t2.*t7.*t355;
t634 = t2.*t7.*t357;
t505 = t504-t634;
t506 = t2.*t81.*(t53-t63);
t507 = t7.*t75.*t357;
t508 = t27.*t43.*t74;
t509 = qd5.*(t507+t508);
t510 = t7.*t44.*t74;
t511 = t7.*t43.*(t53-t63);
t512 = t510+t511;
t513 = qd6.*t512;
t514 = q2d2.*t33.*t357;
t515 = q2d3.*t33.*t357;
t516 = q2d4.*t33.*t357;
t517 = qd2.*t7.*t43.*t81;
t518 = qd3.*t7.*t43.*t81;
t519 = qd4.*t7.*t43.*t81;
t520 = q2d5.*t7.*t43.*t44;
t521 = q2d1.*t7.*t43.*(t68+t69);
t636 = t74.*t81;
t522 = t509+t513+t514+t515+t516+t517+t518+t519+t520+t521-t636;
t523 = t2.*t7.*t43.*(t243-t330);
t524 = t2.*t7.*t44.*t248;
t525 = t2.*t7.*t43.*t295;
t526 = t2.*t7.*t44.*t266;
t527 = t2.*t7.*t44.*t264;
t528 = t2.*t50.*t248;
t529 = t2.*t7.*t44.*t360;
t530 = t2.*t7.*t44.*t344;
t639 = t2.*t7.*t43.*t329;
t531 = t367+t524+t525-t639;
t532 = t2.*t57.*t248;
t533 = t2.*t27.*t44.*t200;
t534 = t2.*t27.*t43.*t334;
t535 = t2.*t7.*t44.*t270;
t536 = t2.*t7.*t43.*t301;
t537 = t2.*t50.*t200;
t538 = t2.*t27.*t266;
t539 = t2.*t27.*t248;
t540 = t2.*t33.*t43.*(3.0./2.5e1);
t541 = t2.*t7.*t43.*t293;
t542 = t2.*t7.*t38.*t248;
t543 = t251+t257+t542-t2.*t70.*t293;
t544 = t2.*t27.*t43.*(1.17e2./1.0e3);
t545 = t2.*t9.*t44.*(5.71e2./1.0e3);
t546 = t2.*t8.*t27.*t43.*(5.71e2./1.0e3);
t547 = t545+t546;
t548 = t2.*t7.*t43.*t314;
t593 = t43.*t293;
t549 = t382-t593;
t550 = t2.*t27.*t270;
t551 = t2.*t27.*t43.*t163;
t552 = t2.*t7.*t43.*t287;
t553 = t2.*t7.*t35.*t248;
t554 = t2.*t7.*t44.*t163;
t555 = t370+t371+t539+t541;
t556 = t2.*t81.*t293;
t557 = t2.*t7.*t43.*t320;
t558 = t2.*t7.*t35.*t200;
t559 = t2.*t9.*t43.*(5.71e2./1.0e3);
t592 = t2.*t8.*t27.*t44.*(5.71e2./1.0e3);
t560 = t559-t592;
t561 = t10.*t27.*t317;
t562 = t44.*t293;
t563 = t388+t562;
t564 = t7.*t10.*t334;
t565 = t10.*t27.*t324;
t566 = t10.*t50.*t293;
t567 = t10.*t33.*t44.*(3.0./2.5e1);
t568 = t7.*t10.*t44.*t293;
t569 = t10.*t57.*t293;
t570 = t10.*t27.*t292;
t571 = t10.*t27.*t200;
t572 = t7.*t10.*t43.*t163;
t573 = t10.*t27.*t329;
t574 = t10.*t27.*t44.*(1.17e2./1.0e3);
t575 = t10.*t50.*(t143+t144-t249-t250);
t576 = t417+t418;
t577 = t44.*t403;
t578 = t27.*t43.*t402;
t579 = t577+t578;
t581 = t43.*t430;
t582 = t44.*t433;
t583 = t581+t582;
t660 = t27.*t43.*t403;
t584 = t420-t660;
t585 = t583.*t584;
t586 = t248.*t256;
t587 = t10.*t293.*t303;
t588 = t586+t587-t10.*t295.*t329;
t589 = t164.*t248;
t590 = t10.*t125.*t293;
t591 = t589+t590-t10.*t295.*t308;
t594 = t428-t2.*t8.*t576;
t595 = t2.*t9.*t576;
t596 = t431+t595;
t597 = t43.*t403;
t674 = t27.*t44.*t402;
t598 = t597-t674;
t599 = t44.*t430;
t675 = t43.*t433;
t600 = t599-t675;
t601 = t27.*t44.*(t381-t421);
t602 = t434+t601;
t603 = t27.*(8.7e1./5.0e2);
t604 = t7.*(t428-t443);
t605 = t603+t604+1.17e2./1.0e3;
t606 = t7.*t10.*t605.*(t381-t421);
t631 = qd3.*(t356-t497);
t632 = qd4.*(t356-t497);
t635 = t495+t496+t506+t631+t632-q2d5.*t505-qd6.*t499-qd5.*t503-t2.*t50.*t74-q2d2.*t2.*t33.*t43.*t44.*2.0-q2d3.*t2.*t33.*t43.*t44.*2.0-q2d4.*t2.*t33.*t43.*t44.*2.0;
t637 = t523+t524+t525+t526;
t638 = t2.*t50.*t266;
t640 = t2.*t57.*t266;
t641 = t2.*t7.*(t381-t421);
t642 = t538+t539+t540+t541;
t643 = t2.*t27.*t344;
t644 = t370+t371+t538+t540;
t661 = t7.*t43.*(3.0./2.5e1);
t645 = t382-t661;
t646 = t2.*t27.*t360;
t647 = t2.*t7.*t81.*(3.0./2.5e1);
t648 = t2.*t7.*t38.*t266;
t649 = t251+t257+t648-t2.*t7.*t70.*(3.0./2.5e1);
t650 = t2.*t602;
t651 = t7.*t44.*(3.0./2.5e1);
t652 = t388+t651;
t653 = t10.*t27.*t44.*t163;
t657 = t10.*t27.*t309;
t654 = t567+t568-t657-t10.*t27.*t295;
t655 = t7.*t10.*t50.*(3.0./2.5e1);
t656 = t10.*t27.*(-t243+t305+t306);
t658 = t7.*t10.*t57.*(3.0./2.5e1);
t659 = t304+t396+t658-t7.*t10.*t38.*t309;
t662 = t243-t330;
t663 = t256.*t266;
t664 = t7.*t10.*t303.*(3.0./2.5e1);
t665 = t663+t664-t10.*t309.*t329;
t666 = t164.*t266;
t667 = t7.*t10.*t125.*(3.0./2.5e1);
t668 = t666+t667-t10.*t308.*t309;
t669 = t248.*t266;
t670 = t10.*t295.*(t243-t330);
t671 = t7.*t10.*t293.*(3.0./2.5e1);
t672 = t669+t670+t671;
t673 = t431+t432;
t676 = t2.*t44.*t57;
t677 = t2.*t43.*t70;
t678 = t2.*t44.*t50;
t679 = t2.*t43.*t81;
t680 = t678+t679;
t681 = t2.*t43.*t256;
t682 = t2.*t44.*t329;
t683 = t681+t682;
t684 = t2.*t43.*t248;
t685 = t684-t2.*t44.*t295;
t686 = t2.*t43.*t266;
t687 = t686-t2.*t44.*t309;
t688 = t2.*t43.*t164;
t689 = t2.*t44.*t70.*(1.17e2./1.0e3);
t690 = t2.*t43.*t57.*(1.17e2./1.0e3);
t691 = t2.*t44.*t81.*(1.17e2./1.0e3);
t692 = t2.*t43.*t50.*(1.17e2./1.0e3);
t693 = t2.*t44.*t125;
t694 = t693-t2.*t7.*t38.*t43.*(1.17e2./1.0e3);
t695 = t544-t2.*t44.*t293;
t696 = t2.*t7.*t35.*t43.*(1.17e2./1.0e3);
t697 = t544-t2.*t44.*t303;
t698 = t544-t2.*t7.*t44.*(3.0./2.5e1);
t699 = t10.*t43.*t303;
t700 = t574+t699;
t701 = t7.*t10.*t43.*(3.0./2.5e1);
t702 = t574+t701;
t703 = t10.*t43.*t125;
t704 = t7.*t10.*t38.*t44.*(1.17e2./1.0e3);
t705 = t10.*t43.*t293;
t706 = t574+t705;
t707 = t43.*t256.*(1.17e2./1.0e3);
t708 = t10.*t44.*t329.*(1.17e2./1.0e3);
t709 = t707+t708;
t710 = t43.*t248.*(1.17e2./1.0e3);
t711 = t710-t10.*t44.*t295.*(1.17e2./1.0e3);
t712 = t43.*t266.*(1.17e2./1.0e3);
t713 = t712-t10.*t44.*t309.*(1.17e2./1.0e3);
t714 = t43.*t164.*(1.17e2./1.0e3);
t715 = t10.*t44.*t308.*(1.17e2./1.0e3);
t716 = t603+t604;
t717 = t53-t63;
Y_l = reshape([qd2.*(t32+t39)+qd3.*(t32+t39)+qd4.*(t32+t39)+qd5.*(t7.*t10.*t38.*t47-t10.*t27.*t30.*t38)-q2d6.*t7.*t10.*t38+q2d1.*t10.*t33.*t38.^2-q2d2.*t7.*t10.*t27.*t38-q2d3.*t7.*t10.*t27.*t38-q2d4.*t7.*t10.*t27.*t38,t453,t453,t453,t10.*t30.*t47,q2d6.*t10-(qd5.*t10.*(t7.*t28.*2.0+t26.*t27.*2.0))./2.0+q2d2.*t10.*t27+q2d3.*t10.*t27+q2d4.*t10.*t27-q2d1.*t7.*t10.*t38-qd2.*t7.*t10.*t35-qd3.*t7.*t10.*t35-qd4.*t7.*t10.*t35,q2d2.*(t66+t67)+q2d3.*(t66+t67)+q2d4.*(t66+t67)+qd6.*(t10.*t30.*t70-t7.*t10.*t38.*t74)-qd2.*t61-qd3.*t61-qd4.*t61+qd5.*(-t10.*t47.*t57+t10.*t27.*t38.*(t53-t63)+t7.*t10.*t30.*t38.*t44+t7.*t10.*t38.*t44.*t75)+q2d6.*t10.*(t56-t64)-q2d5.*t7.*t10.*t38.*t43-q2d1.*t7.*t10.*t38.*t57.*2.0,t466,t466,t466,q2d6.*t10.*t43-t10.*t47.*t54+q2d2.*t10.*t27.*t43+q2d3.*t10.*t27.*t43+q2d4.*t10.*t27.*t43+qd6.*t10.*t30.*t44-qd5.*t10.*t43.*t47+t10.*t30.*t44.*(t29-t41)-q2d1.*t7.*t10.*t38.*t43-qd2.*t7.*t10.*t35.*t43-qd3.*t7.*t10.*t35.*t43-qd4.*t7.*t10.*t35.*t43,q2d5.*t10.*t43+qd2.*t10.*t50+qd3.*t10.*t50+qd4.*t10.*t50+qd6.*t10.*t74-t10.*t30.*t74+q2d1.*t10.*(t56-t64)-q2d2.*t7.*t10.*t44-q2d3.*t7.*t10.*t44-q2d4.*t7.*t10.*t44-qd5.*t10.*t44.*t75,-qd6.*(t2.*t30.*t57+t2.*t7.*t38.*t54)+qd2.*t83+qd3.*t83+qd4.*t83-qd5.*(t2.*t47.*t70+t2.*t27.*t38.*t74+t2.*t7.*t30.*t38.*t43+t2.*t7.*t38.*t43.*t75)+q2d2.*(t71-t2.*t33.*t38.*t43)+q2d3.*(t71-t2.*t33.*t38.*t43)+q2d4.*(t71-t2.*t33.*t38.*t43)+q2d6.*t2.*(t68+t69)-q2d5.*t2.*t7.*t38.*t44-q2d1.*t2.*t7.*t38.*t70.*2.0,t483,t483,t483,q2d6.*t2.*t44+t2.*t47.*t74+q2d2.*t2.*t27.*t44+q2d3.*t2.*t27.*t44+q2d4.*t2.*t27.*t44-qd6.*t2.*t30.*t43-qd5.*t2.*t44.*t47-t2.*t30.*t43.*t75-q2d1.*t2.*t7.*t38.*t44-qd2.*t2.*t7.*t35.*t44-qd3.*t2.*t7.*t35.*t44-qd4.*t2.*t7.*t35.*t44,q2d1.*t2.*(t68+t69)+q2d5.*t2.*t44+qd2.*t2.*t81+qd3.*t2.*t81+qd4.*t2.*t81-t2.*t30.*t54+qd6.*t2.*(t53-t63)+q2d2.*t2.*t7.*t43+q2d3.*t2.*t7.*t43+q2d4.*t2.*t7.*t43+qd5.*t2.*t43.*(t29-t41),-qd5.*(t10.*t44.*t57.*t75+t7.*t10.*t38.*t44.*t54)+qd2.*(t87-t10.*t54.*t60)+qd3.*(t87-t10.*t89.*(t53-t63))+qd4.*(t87-t10.*t89.*(t53-t63))+qd6.*(t10.*t57.*t74-t10.*(t68+t69).*(t53-t63))+q2d1.*t10.*t90.^2+q2d5.*t10.*t43.*(t56-t64)-q2d2.*t7.*t10.*t44.*t57-q2d3.*t7.*t10.*t44.*t57-q2d4.*t7.*t10.*t44.*t57,t494,t494,t494,-qd6.*(t10.*t44.*t54-t10.*t43.*t74)+q2d5.*t10.*t357+qd2.*t10.*t43.*t50+qd3.*t10.*t43.*t50+qd4.*t10.*t43.*t50-t10.*t44.*t54.*t75+q2d1.*t10.*t43.*(t56-t64)-q2d2.*t7.*t10.*t43.*t44-q2d3.*t7.*t10.*t43.*t44-q2d4.*t7.*t10.*t43.*t44-qd5.*t10.*t43.*t44.*t75,t10.*t74.*(t53-t63),q2d5.*(t676+t677)+q2d2.*t92+q2d3.*t92+q2d4.*t92+qd2.*(t94+t95+t96-t2.*t77.*(t53-t63))+qd3.*(t94+t95+t96-t2.*t77.*(t53-t63))+qd4.*(t94+t95+t96-t2.*t77.*(t53-t63))+qd6.*(t2.*t54.*t57.*2.0+t2.*t70.*t74.*2.0)+qd5.*(t2.*t43.*t57.*t75-t2.*t44.*(t68+t69).*(t29-t41)+t2.*t7.*t38.*t43.*t54+t2.*t7.*t38.*t44.*t74)+q2d1.*t2.*(t68+t69).*(t56-t64).*2.0,t495+t496+t506-q2d5.*t505-qd6.*t499-qd5.*t503+qd3.*(t356-t2.*t7.*t44.*t81)+qd4.*(t356-t2.*t7.*t44.*t81)-t2.*t50.*t74-q2d2.*t2.*t33.*t43.*t44.*2.0-q2d3.*t2.*t33.*t43.*t44.*2.0-q2d4.*t2.*t33.*t43.*t44.*2.0,t635,t635,q2d1.*(t676+t677)-q2d2.*t505-q2d3.*t505-q2d4.*t505+qd2.*t680+qd3.*t680+qd4.*t680-qd5.*(t2.*t75.*t355-t2.*t75.*t357)+qd6.*(t2.*t44.*t74.*2.0+t2.*t43.*(t53-t63).*2.0)+q2d5.*t2.*t43.*t44.*2.0+t2.*t44.*t74.*(t29-t41)+t2.*t43.*(t29-t41).*(t53-t63),-t2.*t74.^2+t2.*t717.^2,qd5.*(t43.*(t68+t69).*(t29-t41)-t7.*t38.*t43.*t74)+qd2.*t99+qd3.*t99+qd4.*t99+q2d1.*t100.^2+qd6.*(t54.*t70-t57.*t74)+q2d5.*t44.*(t68+t69)+q2d2.*t7.*t43.*(t68+t69)+q2d3.*t7.*t43.*(t68+t69)+q2d4.*t7.*t43.*(t68+t69),t522,t522,t522,q2d5.*t355+qd6.*(t44.*t54-t43.*t74)+q2d1.*t44.*(t68+t69)+qd2.*t44.*t81+qd3.*t44.*t81+qd4.*t44.*t81-t43.*t74.*t75+q2d2.*t7.*t43.*t44+q2d3.*t7.*t43.*t44+q2d4.*t7.*t43.*t44+qd5.*t43.*t44.*(t29-t41),-t54.*t74,-q2d3.*(t242+t532-t2.*t70.*t295-t2.*t7.*t43.*t153)-q2d4.*(t242+t640-t2.*t70.*t309-t2.*t7.*t43.*t153)-q2d2.*(t242+t366+t2.*t70.*(t305+t306-t7.*t43.*(1.17e2./1.0e3))-t2.*t7.*t43.*t153)+qd4.*(t221+t224+t225+t2.*t74.*(t216+t312-t43.*t114-t27.*t31.*t43.*(1.17e2./1.0e3))-t2.*t89.*t200+t2.*t70.*t317+t2.*t57.*(t203+t276+t277-t27.*t35.*t44.*(1.17e2./1.0e3))-t2.*t275.*(t53-t63))+qd2.*(t221+t224+t225+t2.*t74.*(t216+t322-t43.*t114-t27.*t31.*t43.*(1.17e2./1.0e3))-t2.*t89.*t200+t2.*t70.*t324+t2.*t57.*(t203+t283+t285-t27.*t35.*t44.*(1.17e2./1.0e3))-t2.*t338.*(t53-t63))+q2d5.*(t688+t689+t690+t2.*t44.*t153)+qd3.*(t221+t224+t225-t2.*t89.*t200+t2.*t57.*t264+t2.*t70.*t292+t2.*t74.*t290-t2.*t261.*(t53-t63))+q2d1.*(t2.*t153.*(t68+t69).*2.0+t2.*t57.*t164.*2.0)-qd5.*(-t2.*t54.*t268+t2.*t57.*t270-t2.*t70.*t301+t2.*t74.*t298-t2.*t43.*t75.*t153+t2.*t44.*t75.*t164+t2.*t7.*t38.*t43.*t197+t2.*t7.*t38.*t44.*t200),t537-q2d1.*(t242+t366+t2.*t70.*t329-t2.*t7.*t43.*t308)-qd4.*(t361+t362+t530+t2.*t54.*t373+t2.*t74.*t394-t2.*t7.*t43.*t317)+q2d3.*t531-q2d5.*t683-qd3.*(t361+t362+t527+t2.*t54.*t386+t2.*t74.*(t387+t27.*t44.*t365)-t2.*t7.*t43.*t292)+q2d2.*(t2.*t7.*t44.*t256.*2.0-t2.*t7.*t43.*t329.*2.0)+qd5.*(t533+t534+t535+t536+t2.*t54.*t383+t2.*t74.*t390-t2.*t43.*(t29-t41).*(-t243+t305+t306)+t2.*t44.*t75.*t256)+q2d4.*(t367+t523+t526-t2.*t7.*t43.*t329)-qd2.*(t361+t362+t529-t2.*t7.*t43.*t324)-t2.*t74.*t324+t2.*t7.*t403-t2.*t81.*t334+t2.*t360.*(t53-t63),t537+t641-q2d1.*(t242+t532-t2.*t70.*t295-t2.*t7.*t43.*t308)-qd3.*(t527+t528-t2.*t81.*t295-t2.*t7.*t43.*t292)-qd2.*(t528+t529-t2.*t81.*t295-t2.*t7.*t43.*t324)+q2d2.*t531+q2d4.*t637-q2d5.*t685-qd4.*(t528+t530-t2.*t81.*t295+t2.*t54.*t547+t2.*t74.*t560-t2.*t7.*t43.*t317)+q2d3.*(t2.*t7.*t44.*t248.*2.0+t2.*t7.*t43.*t295.*2.0)+qd5.*(t533+t534+t535+t536+t2.*t54.*t549+t2.*t74.*t563+t2.*t44.*t75.*t248+t2.*t43.*t75.*t295)-t2.*t74.*t292-t2.*t81.*t334+t2.*t264.*(t53-t63),t537+t641-q2d1.*(t242+t640-t2.*t70.*t309-t2.*t7.*t43.*t308)-qd3.*(t527+t638-t2.*t81.*t309-t2.*t7.*t43.*t292)-qd4.*(t530+t638-t2.*t81.*t309-t2.*t7.*t43.*t317)-qd2.*(t529+t638-t2.*t81.*t309-t2.*t7.*t43.*t324)+qd5.*(t533+t534+t2.*t74.*t652+t2.*t645.*(t53-t63)+t2.*t7.*t43.*(t299+t300)+t2.*t7.*t44.*(t269-t43.*(t143+t144-t249-t250))+t2.*t44.*t266.*(t29-t41)+t2.*t43.*(t29-t41).*(t243-t330))+q2d3.*t637-q2d5.*t687+q2d2.*(t367+t523+t526-t639)+q2d4.*(t2.*t7.*t44.*t266.*2.0+t2.*t7.*t43.*t309.*2.0)-t2.*t74.*t317-t2.*t81.*t334+t2.*t344.*(t53-t63),-q2d2.*t683-q2d3.*t685-q2d4.*t687+q2d1.*(t688+t689+t690+t2.*t44.*t308)+qd5.*(t2.*t44.*(t299+t300)-t2.*t43.*t270)+qd3.*(t691+t692+t2.*t43.*t264+t2.*t44.*t292)+qd4.*(t691+t692+t2.*t44.*t317+t2.*t43.*t344)+qd2.*(t691+t692+t2.*t44.*t324+t2.*t43.*t360)+q2d5.*(t2.*t355.*(1.17e2./5.0e2)+t2.*t357.*(1.17e2./5.0e2))-t2.*t54.*t270-t2.*t74.*t301-t2.*t27.*t402-t2.*t44.*t75.*t200-t2.*t43.*t75.*t334,0.0,-q2d1.*(t2.*t70.*t125.*2.0-t2.*t7.*t38.*t164.*2.0)-q2d2.*t369-q2d3.*t543-q2d4.*t649-q2d5.*t694-qd5.*(-t2.*t47.*t164+t2.*t74.*t145+t2.*t70.*t165-t2.*t30.*t268+t2.*t43.*t75.*t125+t2.*t27.*t38.*t200+t2.*t7.*t38.*t270+t2.*t7.*t38.*t43.*t163)+qd3.*(t279+t281+t282-t2.*t81.*(t123+t124-t244)-t2.*t30.*t261+t2.*t70.*t287-t2.*t74.*t288+t2.*t7.*t38.*t264)+qd4.*(t2.*t77.*(t143+t144-t249-t250)+t2.*t314.*(t68+t69)-t2.*t81.*t125-t2.*t30.*t275-t2.*t74.*t310+t2.*t7.*t38.*t344+t2.*t7.*t200.*(t22-t40)+t2.*t7.*t35.*(t146+t147-t201-t202))+qd2.*(t279+t281+t282-t2.*t81.*(t123+t124-t244)-t2.*t30.*t338+t2.*t70.*t320-t2.*t74.*t321+t2.*t7.*t38.*(t203-t278+t283+t285))-qd6.*(t2.*t54.*t125+t2.*t30.*t153+t2.*t57.*t163-t2.*t7.*t38.*t197)-q2d6.*t2.*t164,t558+qd2.*(t374+t557-t2.*t27.*t360-t2.*t7.*t35.*t256)+q2d2.*(t2.*t27.*t256.*2.0+t2.*t7.*t43.*t303.*2.0)-q2d1.*t369+q2d3.*t555+q2d4.*t644-q2d5.*t697+t2.*t436+qd6.*(t554+t2.*t54.*t303+t2.*t30.*t329-t2.*t27.*t334)+qd5.*(t550+t551-t2.*t7.*t200-t2.*t47.*t256+t2.*t74.*t234+t2.*t30.*t383+t2.*t43.*t75.*t303-t2.*t7.*t43.*(t161+t162-t194-t195))+qd4.*(t374+t548-t2.*t27.*t344-t2.*t30.*t373-t2.*t7.*t35.*t256+t2.*t7.*t74.*t229)+qd3.*(t374+t552-t2.*t27.*t264-t2.*t30.*t386-t2.*t7.*t35.*t256+t2.*t7.*t74.*t365)+q2d6.*t2.*t256-t2.*t81.*t163+t2.*t30.*t360-t2.*t74.*t320,t558+t650+q2d3.*(t2.*t27.*t248.*2.0+t2.*t7.*t43.*t293.*2.0)-q2d1.*t543+q2d2.*t555+q2d4.*t642-q2d5.*t695-qd4.*(-t548+t553+t643-t2.*t81.*t293+t2.*t30.*t547+t7.*t8.*t10.*t74.*(5.71e2./1.0e3))+qd3.*(t552-t553+t556-t2.*t27.*t264)-qd2.*(t553-t556-t557+t646)+qd6.*(t554-t2.*t30.*t295+t2.*t54.*t293-t2.*t27.*t334)+qd5.*(t550+t551-t2.*t7.*t200-t2.*t47.*t248+t2.*t74.*t241+t2.*t30.*t549+t2.*t43.*t75.*t293-t2.*t7.*t43.*(t161+t162-t194-t195))+q2d6.*t2.*t248-t2.*t81.*t163+t2.*t30.*t264-t2.*t74.*t287,t558+t650+qd3.*(t552+t647-t2.*t27.*t264-t2.*t7.*t35.*t266)+qd6.*(t554+t2.*t7.*t54.*(3.0./2.5e1)-t2.*t27.*t334-t2.*t30.*(t243-t330))+qd4.*(t548-t643+t647-t2.*t7.*t35.*t266)+qd2.*(t557-t646+t647-t2.*t7.*t35.*t266)+q2d3.*t642+q2d2.*t644-q2d1.*t649-q2d5.*t698+q2d4.*(t2.*t33.*t43.*(6.0./2.5e1)+t2.*t27.*t266.*2.0)+qd5.*(t550+t551+t2.*t27.*t74.*(3.0./2.5e1)-t2.*t7.*t200-t2.*t47.*t266+t2.*t30.*t645+t2.*t7.*t43.*t75.*(3.0./2.5e1)-t2.*t7.*t43.*(t161+t162-t194-t195))+q2d6.*t2.*t266-t2.*t81.*t163+t2.*t30.*t344-t2.*t74.*t314,-q2d1.*t694-q2d3.*t695-q2d2.*t697-q2d4.*t698+qd3.*(t696+t2.*t44.*t287)+qd4.*(t696+t2.*t44.*t314)+qd2.*(t696+t2.*t44.*t320)-qd6.*(t2.*t30.*t44.*(1.17e2./1.0e3)+t2.*t43.*t163)+qd5.*(t2.*t43.*t47.*(1.17e2./1.0e3)-t2.*t44.*t165)+t2.*t74.*(t161+t162-t194-t195)-q2d6.*t2.*t43.*(1.17e2./1.0e3)+t2.*t47.*t200-t2.*t30.*t270-t2.*t43.*t75.*t163+t2.*t7.*t44.*t402,t2.*(t577+t578)-q2d5.*t2.*t43.*(1.17e2./1.0e3)-q2d1.*t2.*t164+q2d3.*t2.*t248+q2d2.*t2.*t256+q2d4.*t2.*t266-qd3.*t2.*t264-qd6.*t2.*t334-qd4.*t2.*t344-qd2.*t2.*t360-t2.*t54.*t163+t2.*t30.*t334+qd5.*t2.*(t269-t43.*(t143+t144-t249-t250)),q2d5.*(t703+t704)-q2d4.*t659-qd3.*(t315+t318-t10.*t50.*t125+t10.*t30.*t290+t10.*t54.*t288+t10.*t57.*t287-t7.*t10.*t35.*t153-t7.*t10.*t38.*t292)-qd4.*(t315+t318-t10.*t50.*t125+t10.*t54.*t310+t10.*t57.*t314+t10.*t30.*t347-t7.*t10.*t35.*t308-t7.*t10.*t38.*t317)-qd2.*(t315+t318-t10.*t50.*t125+t10.*t30.*t342+t10.*t54.*t321+t10.*t57.*t320-t7.*t10.*t35.*t308-t7.*t10.*t38.*t324)+qd5.*(-t10.*(t53-t63).*(t115+t119-t149)+t10.*t47.*t153+t10.*t57.*t165+t10.*t30.*t298-t10.*t44.*(t29-t41).*(t123+t124-t244)+t10.*t27.*t38.*t197+t7.*t10.*t38.*t301-t7.*t10.*t38.*t44.*(t143+t144-t249-t250))-q2d3.*(t302+t304+t569-t7.*t10.*t38.*t295)+qd6.*(-t10.*(t68+t69).*(t143+t144-t249-t250)+t10.*t30.*t164+t10.*t74.*t125+t7.*t10.*t38.*t200)+q2d1.*(t10.*(t56-t64).*(t123+t124-t244).*2.0+t7.*t10.*t38.*t308.*2.0)-q2d2.*(t302+t304+t397+t398)-q2d6.*t10.*t308,t420+t575-qd3.*(t570+t10.*t50.*t303-t10.*t30.*t407-t7.*t10.*t44.*t287+t7.*t10.*t35.*t329-t7.*t10.*(t53-t63).*(t228-t385))+q2d2.*(t10.*t27.*t329.*2.0+t7.*t10.*t44.*t303.*2.0)-q2d5.*t700+q2d3.*(t399+t568+t573-t10.*t27.*t295)+q2d4.*(t399+t567+t656-t10.*t27.*t309)+qd5.*(t564+t653+t10.*t54.*t234-t10.*t27.*t301-t10.*t47.*t329-t10.*t30.*t390-t7.*t10.*t44.*t165+t10.*t44.*t75.*t303)-qd6.*(t571+t572+t10.*t30.*t256+t10.*t74.*t303)-qd4.*(t391+t395+t561-t10.*t30.*t394-t7.*t10.*t54.*t229-t7.*t10.*t44.*t314)-qd2.*(t391+t395+t565-t7.*t10.*t44.*t320)-q2d1.*(t304+t396+t397+t398)+q2d6.*t10.*(-t243+t305+t306)+t10.*t30.*t324-t10.*t54.*t320-t27.*t43.*t403-t7.*t10.*t35.*t334,t420+t575+qd5.*(t564-t10.*t27.*t301+t10.*t47.*t295-t10.*t30.*t563+t10.*t241.*(t53-t63)-t7.*t10.*t44.*t165+t10.*t44.*t293.*(t29-t41)+t10.*t27.*t44.*(t143+t144-t249-t250))-q2d3.*(t10.*t27.*t295.*2.0-t7.*t10.*t44.*t293.*2.0)+q2d4.*t654-q2d5.*t706+q2d2.*(t399+t568+t573-t10.*t27.*t295)-qd3.*(t566+t570-t7.*t10.*t35.*t295-t7.*t10.*t44.*t287)-qd2.*(t565+t566-t7.*t10.*t35.*t295-t7.*t10.*t44.*t320)-qd6.*(t571+t572+t10.*t30.*t248+t10.*t74.*t293)-q2d1.*(t304+t396+t569-t7.*t10.*t38.*t295)-qd4.*(t561+t566-t10.*t30.*t560-t7.*t10.*t35.*t295-t7.*t10.*t44.*t314+t2.*t7.*t8.*t10.*t54.*(5.71e2./1.0e3))-q2d6.*t10.*t295+t10.*t30.*t292-t10.*t54.*t287-t27.*t43.*t403-t7.*t10.*t35.*t334,t420+t575-t660+q2d3.*t654-q2d1.*t659-q2d5.*t702+q2d2.*(t399+t567+t656-t657)-qd3.*(t570+t655-t7.*t10.*t44.*t287-t7.*t10.*t35.*t309)-qd4.*(t561+t655-t7.*t10.*t35.*t309-t7.*t10.*t44.*t314)-qd2.*(t565+t655-t7.*t10.*t35.*t309-t7.*t10.*t44.*t320)-qd6.*(t571+t572+t7.*t10.*t74.*(3.0./2.5e1)+t10.*t30.*t266)+q2d4.*(t10.*t33.*t44.*(6.0./2.5e1)-t10.*t27.*t309.*2.0)+qd5.*(t564+t653-t10.*t27.*(t299+t300)+t10.*t27.*t54.*(3.0./2.5e1)+t10.*t47.*t309-t10.*t30.*t652+t7.*t10.*t44.*t75.*(3.0./2.5e1)-t7.*t10.*t44.*(t161+t162-t194-t195))-q2d6.*t10.*t309+t10.*t30.*t317-t10.*t54.*t314-t7.*t10.*t35.*t334,q2d1.*(t703+t704)-qd3.*(t10.*t43.*t287-t7.*t10.*t35.*t44.*(1.17e2./1.0e3))-qd4.*(t10.*t43.*t314-t7.*t10.*t35.*t44.*(1.17e2./1.0e3))-qd2.*(t10.*t43.*t320-t7.*t10.*t35.*t44.*(1.17e2./1.0e3))-q2d2.*t700-q2d4.*t702-q2d3.*t706+qd6.*(t10.*t30.*t43.*(1.17e2./1.0e3)-t10.*t44.*t163)+qd5.*(t10.*t44.*t47.*(1.17e2./1.0e3)+t10.*t43.*t165)+t10.*(t53-t63).*(t161+t162-t194-t195)+t10.*t30.*(t299+t300)-q2d6.*t10.*t44.*(1.17e2./1.0e3)-t10.*t47.*t334-t7.*t43.*t402-t10.*t44.*t75.*t163,-t597+t674+t10.*t74.*(t143+t144-t249-t250)+q2d2.*t10.*(-t243+t305+t306)-q2d5.*t10.*t44.*(1.17e2./1.0e3)-q2d3.*t10.*t295-q2d1.*t10.*t308-q2d4.*t10.*t309-qd6.*t10.*t200-qd3.*t10.*t292-qd5.*t10.*t301-qd4.*t10.*t317-qd2.*t10.*t324+t10.*t30.*t200,q2d5.*(t714+t715)-qd2.*(t200.*t338-t360.*(t146+t147-t201-t202)+t10.*t125.*t320+t10.*t163.*t321-t10.*t308.*t324-t10.*t334.*t342)+qd5.*(t200.*t268-(t269-t43.*(t143+t144-t249-t250)).*(t146+t147-t201-t202)-t10.*t334.*(t296+t297)+t10.*t125.*t165+t10.*t301.*t308-t10.*(t115+t119-t149).*(t143+t144-t249-t250))-q2d2.*t411-q2d3.*t591-q2d4.*t668+q2d1.*(t10.*t308.^2+t10.*t326.^2+t325.^2)+qd3.*(t164.*t264-t200.*t261-t10.*t287.*(t123+t124-t244)+t10.*t292.*t308+t10.*t290.*t334-t10.*(t272+t273-t336).*(t143+t144-t249-t250))-qd4.*(t200.*t275-t164.*t344+t10.*t125.*t314+t10.*t163.*t310-t10.*t308.*t317-t10.*t334.*t347)+qd6.*(t164.*t197-t200.*t308-t10.*t334.*(t146+t147-t201-t202)+t10.*t200.*t308),t585+t606-qd4.*(t200.*t373+t256.*t344-t10.*t303.*t314+t10.*t317.*t329+t10.*t334.*t394-t7.*t10.*t163.*t229)-qd3.*(t256.*t264+t200.*t386-t10.*t287.*t303+t10.*t292.*t329+t10.*t334.*t407-t7.*t10.*t163.*t365)+t579.*(t44.*t440+t27.*t43.*t442)-q2d1.*t411+q2d3.*t588+q2d4.*t665-q2d5.*t709+t200.*t360+q2d2.*(t10.*t303.^2+t10.*t404.^2+t256.^2)+qd5.*(t256.*t270+t200.*t383-t10.*t303.*(t161+t162-t194-t195)-t10.*(t299+t300).*(-t243+t305+t306)+t10.*t163.*t234+t10.*t334.*t390)+qd6.*(t200.*t329-t256.*t334-t10.*t200.*(-t243+t305+t306)+t10.*t256.*t334)-qd2.*(t256.*t360-t10.*t303.*t320+t10.*t324.*t329)+t10.*t598.*(t43.*t440-t27.*t44.*t442)-t10.*t163.*t320-t10.*t324.*t334-t10.*t436.*t600+t10.*t33.*t402.*t442,t585+t606+qd5.*(t200.*t549+t248.*(t269-t43.*(t143+t144-t249-t250))+t10.*t241.*(t143+t144-t249-t250)+t10.*t295.*(t299+t300)-t10.*t165.*t293+t10.*t334.*t563)-t579.*(t44.*t594+t27.*t43.*t596)-qd6.*(t200.*t295+t248.*t334-t10.*t200.*t295-t10.*t248.*t334)+q2d2.*t588-q2d1.*t591+q2d4.*t672-q2d5.*t711+t200.*t264+q2d3.*(t10.*t293.^2+t10.*t295.^2+t248.^2)-qd4.*(t248.*t344+t200.*t547-t10.*t293.*t314-t10.*t295.*t317+t10.*t334.*t560+t2.*t7.*t8.*t10.*t163.*(5.71e2./1.0e3))+qd3.*(-t248.*t264+t10.*t287.*t293+t10.*t292.*t295)+qd2.*(-t248.*t360+t10.*t293.*t320+t10.*t295.*t324)-t10.*t598.*(t43.*t594-t27.*t44.*t596)-t10.*t163.*t287-t10.*t292.*t334-t10.*t600.*t602-t10.*t33.*t402.*t596,t585+t606-t579.*(t44.*t429+t27.*t43.*t673)-qd6.*(t200.*t309+t266.*t334-t10.*t200.*t309-t10.*t266.*t334)+q2d2.*t665-q2d1.*t668+q2d3.*t672-q2d5.*t713+t200.*t344+qd5.*(t200.*t645+t266.*(t269-t43.*(t143+t144-t249-t250))+t10.*t27.*(t143+t144-t249-t250).*(3.0./2.5e1)-t7.*t10.*t165.*(3.0./2.5e1)+t10.*t334.*t652+t10.*(t299+t300).*(t243-t330))+qd3.*(-t264.*t266+t7.*t10.*t287.*(3.0./2.5e1)+t10.*t292.*t309)+qd4.*(-t266.*t344+t7.*t10.*t314.*(3.0./2.5e1)+t10.*t309.*t317)+qd2.*(-t266.*t360+t7.*t10.*t320.*(3.0./2.5e1)+t10.*t309.*t324)+q2d4.*(t10.*t33.*(9.0./6.25e2)+t10.*t662.^2+t266.^2)-t10.*t598.*(t43.*t429-t27.*t44.*t673)-t10.*t163.*t314-t10.*t317.*t334-t10.*t600.*t602-t10.*t33.*t402.*t673,q2d1.*(t714+t715)+qd3.*(t43.*t264.*(1.17e2./1.0e3)+t10.*t44.*t292.*(1.17e2./1.0e3))-qd5.*(t43.*t270.*(1.17e2./1.0e3)-t10.*t44.*t301.*(1.17e2./1.0e3))+qd4.*(t43.*t344.*(1.17e2./1.0e3)+t10.*t44.*t317.*(1.17e2./1.0e3))+qd2.*(t43.*t360.*(1.17e2./1.0e3)+t10.*t44.*t324.*(1.17e2./1.0e3))-qd6.*(t44.*t200.*(1.17e2./1.0e3)-t43.*t334.*(1.17e2./1.0e3)-t10.*t44.*t200.*(1.17e2./1.0e3)+t10.*t43.*t334.*(1.17e2./1.0e3))-q2d2.*t709-q2d3.*t711-q2d4.*t713-t200.*t270+q2d5.*(t357.*1.3689e-2+t10.*t355.*1.3689e-2)+t43.*t716.*(t577+t578)+t10.*(t161+t162-t194-t195).*(t143+t144-t249-t250)-t10.*t301.*t334+t7.*t10.*t402.*t430-t7.*t43.*t402.*t583-t10.*t27.*t402.*t605-t10.*t44.*t598.*t716-t7.*t10.*t44.*t402.*t600,t600.*(t577+t578)+t200.*t334-t583.*t598-t10.*t200.*t334-t10.*t579.*t600+t10.*t583.*(t597-t674)],[6,10]);

author: shamilmamedov


text:
function [ s1, s2 ] = s_swap ( s1, s2 )

%*****************************************************************************80
%
%% S_SWAP swaps two strings.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S1, S2, two strings.
%
%    Output, string S1, S2, the interchanged strings.
%
  s3 = s1;
  s1 = s2;
  s2 = s3;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% function [B, Bb, Bw, pc_b, sc_b, pc_w, sc_w] = mlpcr2(X,Y,varargin)
%
% Multilevel PCR v3.0 designed for prediction on datasets with multiple
% subjects and multiple observations per subject. Unlike version 2 this
% supports full mixed effects modeling. See randSlope and randInt
% arguments.
%
% Identifies within block (subject) eigenvectors with optional balancing
% across subjects. Computes loadings of full dataset on these within
% subject eigenvectors (so they can vary both within and between blocks).
% Computes a separate PCA on the residual and obtains a second set of
% orthogonal components and scores which only vary between blocks. Performs
% regression of between and within loadings (jointly) on pain outcome, and
% projects regression coefficients back to voxel space using within and
% between eigenvectors.
%
% regression is OLS (default) or moore-penrose pseudoinverse (rank
% deficient data), not some mixed effects thing. Within subject PCA is
% rate limiting step.
%
% If the defaults are used then the result is identical to PCR, only you
% get both within and between subject predictive models in addition to the
% full model. Different results may be obtained with hyperparameter
% optimization or concensus PCA enabled.
%
% Input ::
%
%   X           - n x p data matrix
%
%   Y           - n x 1 outcome vector
%
%   'subjIDs    - n x 1 vector (ideal) or cellstr (throws warning, but
%                   fine) indicating group affilitations. Subjects must be
%                   adjacent. Don't intermix blocks with one another. This
%                   function was not designed for this scenario and will
%                   break even with equivalently intermixed block labels.
%
% Optional Input ::
%
%   'numcomponents'
%               - 1 x 2 numeric vector indicating number of PCA dimensions
%                   to retain at the between and within levels
%                   (respectively). Use Inf to autoselect based on
%                   available degrees of freedom. default: [Inf, Inf]
%                 Note 1: bayesopt from the Statistics and Machine Leraning
%                   matlab toolbox provides an elegant way to optimize
%                   these variables.
%                 Note 2: More dimensions result in less biased solutions,
%                   which is a good reason for using the default [Inf,Inf].
%                   However task driven variance is going to be captured by
%                   within components, so [0, Inf] is also a sensible
%                   default if you wish to treat trait differences as
%                   noise.
%
%   'cpca'      - followed by 0/1 to indicate whether or not concensus PCA
%                   (Westerhuis, et al. 1998) should be used in place of
%                   standard pca. If concensus PCA is enabled eigenvectors
%                   will be selected such that variance is explained
%                   equally across all blocks. Otherwise eigenvectors will
%                   best represent blocks with the most observations.
%                 Note 1: for optimization of dimension hyperparameters
%                   specify a custom loss function that also balances the
%                   weight of each block. fmri_data/predict won't do this
%                   and will consequently fight against CPCA. A future
%                   update may fix this, and if so this note should be
%                   removed.
%                 Note 2: In principle concensus PCA could be implemened
%                   for traditional PCR too. It just hasn't been.
%
%   'randInt'   - Fit a random intercept model. Default=false
%
%   'randSlope' - Fit a random slope model (within effects only).
%                   Default=false
%
%   'fitlmeOpts'
%               - Cell array of options to passthrough to fitlme. Run with
%               {'CovariancePattern','isotropic'} by default.
%
% Output ::
%
%   B           - B(1) is intercept, B(2:end) are regression weights in X
%                   space.
%
%   Bb          - Bb(1) is intercept. Bb(2:end) are regression weights of
%                   between block variance components in X space.
%
%   Bw          - Bw(1) is zero. Bb(2:end) are regression weights of within
%                   block variance components in X space.
%
%   pc_b        - between eigenvectors
%
%   sc_b        - scores on between eigenvectors
%
%   pc_w        - within eigenvectors
%
%   sc_w        - scores on within eigenvectors
%
%
% Version History ::
%
%   MLPCR was originally developed using mixed effects models (version 1,
%   Petre, et al., 2019), however model fitting was prohibitively slow and
%   consequently performance characteristics could not be fully evaluated.
%   Version 2.0 substitutes (optionally weighted) OLS for much faster
%   convergence, and minimaly different weight fitting. Version 1 is much
%   more flexible than version 2, and is still available in CanlabPrivate
%   for potential future development. It has been removed from CanlabCore.
%   Version 3.0 reintroduces mixed effects modelling as a non-default
%   option, but abides by better coding practices, following in the
%   footsteps of version 2.0.
%
% References ::
%
%   Petre B, Woo W, Losin E, Eisenbarth H, Wager TD. (2019) Separate within
%       -subject and individual-difference predictions with multilevel
%       MVPA. Society for Neuroscience, San Diego, CA. (included with
%       canlabCore mlpcr library as pdf).
%
%   Westerhuis J, Kourti T, MacGregor J. (1998) Analysis of Multiblock and
%       Hierarchical PCA and PLS Methods. Journal of Chemometrics 12(5).
%
%
% Designed and writen by Bogdan, 5/4/2020
%
%
%
% ToDo:
% - Enable between or within dimension retention priority for bootstrapping
%   (allowing the user to force retention of one or the other even if
%   eigenvalue rank doesn't justify it)
% - passthrough options to higher level cv_mlpcr, cv_mlpcr_bt and
%   cv_mlpcr_wi scripts. Concensus PCA should also use concensus cv_err,
%   cv_mlpcr_wi and cv_mlpcr_bt should have within and between priority
%   (respectively) by default.

function [B, Bb, Bw, pc_b, sc_b, pc_w, sc_w, b] = mlpcr3(X,Y,varargin)
    subjIDs = [];
    wiDim = Inf;
    btDim = Inf;
    cpca = 0;
    randInt = 0;
    randSlope = 0;
    fitlmeOpts = {};
    for i = 1:length(varargin)
        if ischar(varargin{i})
            switch varargin{i}
                case 'subjIDs'
                    subjIDs = varargin{i+1}(:);
                case 'numcomponents'
                    nc = varargin{i+1};
                    btDim = nc(1);
                    wiDim = nc(2);
                case 'cpca'
                    cpca = varargin{i+1};
                case 'randInt'
                    randInt = varargin{i+1};
                case 'randSlope'
                    randSlope = varargin{i+1};
                case 'fitlmeOpts'
                    fitlmeOpts = varargin{i+1};
            end
        end
    end

    if randInt || randSlope
        if ~ismember('CovariancePattern',fitlmeOpts)
            fitlmeOpts = [fitlmeOpts, {'CovariancePattern', 'isotropic'}];
        end
    end

    % check for intercept and remove if present
    %{
    intercept = ones(1,size(X,2))';
    wh = cellfun(@(x) isequal(x,intercept),num2cell(X,1));
    if ~isempty(wh)
        warning('Your input X should not contain an intercept. Removing it. Model parameters will corespond to intercept free design')
        X(:,wh) = [];
    end
    %}

    % we need adjacent subjIDs, so let's ensure that
    [subjIDs, newOrder] = sortrows(subjIDs(:));
    [~,origOrder] = sort(newOrder);

    X = X(newOrder,:);
    Y = Y(newOrder);

    if isempty(subjIDs)
        error('Cannot perform multilevel PCR without a subjIDs identifier');
    end

    [~, grp_exemplar, subjIDs] = unique(subjIDs,'rows','stable');
    uniq_grp = unique(subjIDs);

    % get centering and expansion matrices
    n_grp = length(uniq_grp);
    cmat = cell(1,n_grp);
    emat = cell(1,n_grp);
    sf = []; % scale factor for imbalanced datasets
    for i = 1:n_grp
        this_grp = uniq_grp(i);
        this_n = sum(this_grp == subjIDs);
        cmat{i} = eye(this_n) - 1/this_n;
        emat{i} = ones(this_n,1);
        sf = [sf(:); 1/sqrt(this_n)*ones(this_n,1)];
    end
    cmat = blkdiag(cmat{:});
    emat = blkdiag(emat{:});
    if ~cpca
        sf = ones(size(sf));
    end

    % compute preliminary within fractions
    Xw = cmat*X;

    if wiDim > 0
        % determine within dimension retention
        if wiDim > length(subjIDs) - length(uniq_grp)
            if wiDim < Inf % something user supplied was too big
                warning('Max wiDim exceeds max df, reseting wiDim to %d',length(subjIDs) - length(uniq_grp));
            end

            wiDim = min(length(subjIDs) - length(uniq_grp), size(X,2));
        end

        % Get concensus PCA solution, which weighs each block equally
        % requires Matlab 2016b or later perform operation across all
        %   columns like this with sf
        [pc_w,~,~] = svd((sf.*Xw)', 'econ');
        pc_w = pc_w(:,1:wiDim);

        sc_w = X*pc_w;

        % note: sc_w*pc_w' ~= Xw, it pulls scores from the entire dataset,
        % 	Some of the between fraction may vary along the within
        % 	components, and we want to pull that into the within scores.
        % 	sc_w*pc_w' is not necessarily mean zero within subject.

        % modified Xb, invariant with subject, but not exactly the subject
        %   ean either due to the missing variance along pc_w
        Xr = X - sc_w*pc_w';
        Xb = Xr - cmat*Xr; % if we don't have full wiDim then there's residual within variance to remove
        Xb = Xb(grp_exemplar,:);
    else
        [pc_w, sc_w] = deal([]);

        Xb = X - Xw;
        Xb = Xb(grp_exemplar,:);
    end

    if btDim > 0
        % get between components from residual between fraction
        if btDim > length(uniq_grp) - 1
            if btDim < Inf % something user supplied was too big
                warning('Max btDim exceeds max df, reseting btDim to %d',length(uniq_grp) - 1);
            end

            btDim = min(length(uniq_grp) - 1, size(Xb,2));
        end

        [pc_b,~,~] = svd(scale(Xb,1)','econ');
        pc_b = pc_b(:,1:btDim);
        sc_b = Xb*pc_b;
    else
        [pc_b, sc_b] = deal([]);
    end

    if ~isempty(sc_b)
        sc_b = emat*sc_b;
        bDim = 1:size(sc_b,2);
        wDim = (1:size(sc_w,2)) + bDim(end);
    else
        bDim = [];
        wDim = 1:size(sc_w,2);
    end
    sc = [sc_b, sc_w];
    pc = [pc_b, pc_w];

    % 3/8/13 TW solution from cv_pcr code to use numcomps, because sc is
    %   not always full rank during bootstrapping.
    % Modified by Bogdan for use in mlpcr to keep componens with the
    %   highest ranking eigenvectors (not necessarily sequential in mlpcr
    %   because of within/between stratification)
    % ToDo: create a bootstrap priority option to allow the user to specify
    %   whether they want to prefer between or within component retention
    %   regardless of eigenvalue size.
    if rank(sc) == size(sc,2)
        numcomps = rank(sc);
        retainComps = 1:numcomps;

        bDim(~ismember(bDim,retainComps)) = [];
        wDim(~ismember(wDim,retainComps)) = [];
    elseif rank(sc) < size(sc,2)
        error('You will run into problems when reconstructing Bw at the end. Fix the wDim and bDim indexing');
        numcomps = rank(sc)-1;
        [~,compRank] = sort(var(sc),'descend');
        retainComps = sort(compRank(1:numcomps));

        [bDimnew, wDimnew] = deal(zeros(length(compRank),1));

        bDim(~ismember(bDim,retainComps)) = [];
        wDim(~ismember(wDim,retainComps)) = [];

        if ~any(ismember(bDim,retainComps)) && btDim > 0
            warning('All between dimensions dropped due to rank deficiency');
        end
        if ~any(ismember(wDim,retainComps)) && wiDim > 0
            warning('All within dimensions dropped due to rank deficiency');
        end

        sc_w = sc(:,wDim);
        sc_b = sc(:,bDim);
        pc_w = pc(:,wDim);
        pc_b = pc(:,bDim);

        bDimnew(bDim) = 1;
        bDimnew = bDimnew(retainComps);
        bDim = find(bDimnew);

        wDimnew(wDim) = 1;
        wDimnew = wDimnew(retainComps);
        wDim = find(wDimnew);
    end

    xx = [ones(size(Y, 1), 1) sc(:, retainComps)];

    if ~randInt && ~randSlope % fit a LS or weighted LS model
        if rank(xx) <= size(sc, 2)
            % compute (optional: weighted) pseudoinverse if not full rank
            tol = max(size(xx))*eps(norm(xx));
            [u,s,v] = svd(sf.*xx,'econ');
            s = diag(s);
            s(s>tol) = 1./s(s>tol);
            s = diag(s);
            pinv_xx = v*s*u';

            b = pinv_xx * (sf.*Y);
        else
            b = inv(xx'*diag(sf.^2)*xx)*xx'*diag(sf.^2)*Y;
        end
    else % fit random effects model
        if rank(xx) <= size(sc,2)
            warning(['Rank deficient data found. This is not supported with mixed effects models.',...
                'If you''re bootstrapping consider using the fixed effects approach instead']);
        end
        xxRE = [];
        if randInt
            xxRE = [xxRE, ones(size(Y,1),1)];
        end
        if randSlope
            xxRE = [xxRE, sc_w];
        end

        %m = fitlmematrix(double(sf.*xx), sf.*Y, double(sf.*xxRE), categorical(subjIDs), fitlmeOpts{:});
        %try
            m = fitlmematrix(double(xx), Y, double(xxRE), categorical(subjIDs), fitlmeOpts{:});
        %catch
        %    keyboard
        %end
        b = m.fixedEffects;
    end

    if ~isempty(retainComps)
        B = [b(1); pc(:,retainComps)*b(2:end)];
    else
        B = [b;  zeros(size(X,2),1)];
    end

    if isempty(bDim)
        Bb = [b(1); zeros(size(X,2),1)];
    else
        %Bb = [b(1); pc_b(:,bDim)*b(bDim + 1)];
        Bb = [b(1); pc(:,bDim)*b(bDim + 1)];
    end

    if isempty(wDim)
        Bw = [0; zeros(size(X,2),1)];
    else
        %Bw = [0; pc_w(:,wDim)*b(wDim + 1)];
        Bw = [0; pc(:,wDim)*b(wDim+1)];
    end

    if ~isempty(sc_b), sc_b = sc_b(origOrder,:); end
    if ~isempty(sc_w), sc_w = sc_w(origOrder,:); end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore